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PREFACE 


This first Annual Report covers the initial effort on the NASA 
HOST program titled "Constitutive Modeling for Isotropic Materials" 
conducted under Contract NAS3-23925. The NASA program manager for 
this project is Mr. Albert Kaufman. The program manager at Southwest 
Research Institute is Dr. Ulric Lindholm. Contributors to this report 
are Dr. Kwai Chan and Mr. Andrew Nagy of SwRI, Messrs. Jeff Hill and 
R. M. Weber of Pratt & Whitney Aircraft, and Dr. Kevin Walker and Prof. 
Sol R. Bodner, consultants. 
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1.0 INTRODUCTION 


The objective of the present program is to develop a unified con- 
stitutive model for finite-element structural analysis of turbine engine 
hot section components. This effort constitutes a different approach for 
non-linear finite-element computer codes which have heretofore been based 
on classical inelastic methods. The unified constitutive theory to be 
developed will avoid the simplifying assumptions of classical theory and 
should more accurately represent the behavior of superalloy materials 
under cyclic loading conditions and high temperature environments. This 
class of constitutive theory is characterized by the use of kinetic 
equations and internal variables with appropriate evolutionary equations 
for treating all aspects of inelastic deformation including plasticity, 
creep and stress relaxation. Model development is directed toward iso- 
tropic, cast nickel-base alloys used for air-cooled turbine blades and 
vanes. Recent studies have shown that this approach is particularly 
suited for determining the cyclic behavior of superalloy type blade and 
vane materials and is entirely compatible with three dimensional inelastic 
finite-element formulations. More efficient and accurate inelastic 
analysis of hot section components--turbi ne blades, turbine vanes, com- 
bustor liners and seals-- fabricated from "age hardenable" isotropic super- 
alloy materials will be realized as the result of these developments. 

The program is being conducted in two phases. A basic program 
(Tasks A through I) and an optional follow-on program (Tasks J through M). 

In the Basic Program of twenty six months' duration, a unified consti- 
tutive model will be developed for the prediction of the structural re- 
sponse of isotropic materials for temperatures and strain range character- 
istics of cooled turbine vanes in advanced gas turbine engines (Task A). 

A data base of uniaxial and multiaxial material properties required for 
the constitutive model development will be obtained for the base material 
(Tasks C and E). The constitutive model will then be incorporated into 
a finite-element computer code (Task D). An evaluation will be made of 
the capability of the analytical method to predict the structural response 
for multiaxial stress states (Task E) and nonisothermal conditions by con- 
ducting thermomechanical loading and benchmark notch verification experi- 
ments and analysis (Task F). As a final evaluation of the analytical model, 
a structural analysis will be performed for a hot section component fabri- 
cated from the base material for simulated engine operating conditions 
(Task G). In the optional program material property test procedures will 
be further developed to minimize the amount of testing required, and to 
study the possibility for estimating the material model constants from 
conventional property data (Task J). Further development of the model 
will be undertaken to consider thermal history effects and to correct any 
deficiencies indicated in the model or the computational algorithms in 
the code (Task K) . In addition, the constitutive model development will 
be verified for an alternate material (Task L). 
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The work under this program is being conducted as a joint effort 
between Pratt & Whitney Aircraft (PWA) and Southwest Research Institute 
(SwRI) with technical assistance from Prof. Sol R. Bodner and Dr. Kevin 
Walker in the area of constitutive model development. The work and data 
base generated is being coordinated closely with another NASA funded 
HOST program at PWA (NAS3-23288) to develop advanced life prediction 
techniques for isotropic superalloy vane and blade materials. 

Progress to date on this program has included completion of the 
review of unified constitutive models (Task A), substantial effort on 
specimen fabrication (Task B), uniaxial testing (Task C), initiation of 
effort on the implementation of models in MARC finite-element code (Task D), 
and multiaxial testing (Task E). A report on technical progress and future 
plans is given in the following sections. 
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2.0 TASK A. REVIEW AND SCREENING OF CANDIDATE CONSTITUTIVE MODELS 

2.1 Literature Survey 

A literature survey has been conducted to assess the state-of-the-art 
of time-temperature dependent elastic-viscoplastic constitutive theories 
which are based on the unified approach. This class of constitutive theories 
is characterized by the use of kinetic equations and internal variables with 
appropriate evolutionary equations for treating all aspects of inelastic de- 
formation including plasticity, creep, and stress relaxation. The review 
identifies more than ten such unified theories which are shown to satisfy the 
uniqueness and stability criteria imposed by Drucker's postulate and Ponter's 
inequalities. These theories are compared on the basis of the types of flow 
law, kinetic equation, evolutionary equation of the internal variables, and 
treatment of temperature dependence. The similarities and differences of 
these theories are first outlined in terms of mathematical formulations and 
then illustrated by comparisons of theoretical calculations with experimental 
results which include monotonic stress-strain curves, cyclic hysteresis loops, 
creep and stress relaxation rates, as well as thermomechanical loops. Numer- 
ical methods used for integrating these stiff time-temperature dependent con- 
stitutive equations are also reviewed. 

Because of its length, this review has been included as Appendix A of 
this report. 

2.2 Model Selection for Subsequent Study 

As a result of the literature survey and based upon SwRI recommendations, 
two constitutive models have been selected and approved for further study. 

These are the models of Bodner and Partom [1] and of Walker [2]. These two 
models are both representati ve of the class of unified models considered in 
the review process but differ significantly in the choice of particular 
functional forms for the basic flow law, the kinetic relationship, the para- 
meter used as a measure of hardening and the evolutionary equations for the 
internal variables describing strain hardeninq. Thus, a direct comparison 
between these two models and the experiments should illustrate well the conse- 
quences of a wide range in constitutive modeling approach. It is also signi- 
ficant that both models have already found significant application to analysis 
of gas turbine material s and to hot section components. Therefore, they are 
further along in their development and evaluation than most of the other 
comparable models. 

2. 3 Identification of Model Modifications 

. As . a result .of the literature survey and of direct experience in 
working with existing models, three specific areas have been identified which 
will potentially require modification to the models as they currently exist 
These areas are listed and discussed below. 
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1. Mul tiaxial , non-proportional load or strain histories . There is 
existing data which indicates that the rate of strain hardening varies 
with a change in direction of the strain rate vector; i.e., when the ratio 
of the plastic strain increments is not held constant as in a uniaxial or 
proportional loading history. This has been demonstrated for both monotoni- 
cally increasing and cyclic loading. This poses potential problems in 
formulation of the evolutionary equations for the hardening variables which 
will require some measure of the rate of change in direction of the strain 
increments and the functional inclusion of this measure in the hardening 
relations. Additional practical difficulties will accrue if the material 
constants associated with the non-proportional hardening cannot be obtained 
from uniaxial tests, since the availability of multiaxial test facilities is 
low and the generation of such data is quite expensive relative to uniaxial 
testing. 

2. Non-isothermal strain cycles . It is not evident at this time that 
isothermal test data will be sufficient input for a model adequate to handle 
non-isothermal problems. At elevated temperatures, metallurgical processes 
affecting strength may be occurring which depend not only on the current 
temperature but upon time at temperature or rate of change of temperature. 
Such effects are generally not incorporated in current models. Again, a 
requirement for non-isothermal testing in order to determine constitutive 
model constants would be a practical limitation. 

3. Description of hardening over a wide strain range . Some current 
experience indicates a difficulty in developing hardening behavior which is 
accurate for both small (near yield) and large (saturation) plastic strain 
values. Multiple functions may be required for a general model. The corres- 
ponding increase in number of material constants required is the penalty 
incurred. 

2.4 References 

1. S. R. Bodner and Y. Partom, ASME J. of Applied Mechanics, Vol . 42, 

p. 385, 1975. 

2. K. P. Walker, NASA Contract Report NASA CR 165533, 1981. 
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3.0 TASK B. SPECIMEN FABRICATION 


3.1 Material Selection 


A B1900+Hf was selected as the model material for study in the base 
program and MAR-M247 the material for consideration in the optional tasks. 
Both materials are extensively utilized in the gas turbine industry and 
considerable benefit will derive from their characterization. A major 
factor in the selection of B1900+Hf is the availability of cyclic and 
monotonic constitutive data from a concurrent NASA HOST program, NAS3- 23288, 
conducted by PWA to develop life prediction analysis methodology for iso- 
tropic blades and vanes. Selection of the same material, material process- 
ing, and specimen configurations in both efforts significantly enhances the 
data base available for life prediction and constitutive model development. 

Sufficent material of each type was procured to satisfy specimen 
requirements in the program. 

The B1900+Hf material was part of a single heat, designated W-0098, 
obtained from Certified Alloy Products Inc., Long Beach, California. The 
chemical composition of this heat is compared to nominal specification 
in Table 3.1 . 


TABLE 3.1 

CHEMICAL COMPOSITION OF B1900+Hf (HEAT W-0098) 


El ement 

Nominal (%) 

Heat W-0098 

C 

0.11 

0.09 

Cr 

8.0 

7.72 

Co 

10.0 

9.91 

Mo 

6.0 

5.97 

A1 

6.0 

6.07 

Ta 

4.3 

4.21 

Ti 

1 .0 

0.99 

B 

0.015 

0.016 

Zr 

0.08 

0.04 

Fe 

0.35* 

0.17 

W 

0.1* 

0.04 

Cb 

0.1* 

0.08 

Bi 

0.5 ppm 

0.1 

Pb 

10.0 ppm 

0.1 

Hf 

1.5 

1.19 

Ni 

Remai nder 

Remai nder 


*Maximum 
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Casting configurations, pour and mold temperatures for the B1900+Hf 
specimens were selected to assure that grain size, material structure and 
integrity match, as closely as possible, those in the PWA test specimens. 
Goals for this phase of the fabrication effort were: 

1. A grain size of ASTM No. 1 to 2 in the Gage section when 
measured using standard procedures. 

2. A porosity-free casting, see Figure 3.1. 

3. A y‘ size of .9 ym in the fully heat treated condition. 

Following casting, all B1900+Hf test bars were fully heat treated 
according to the following schedule: 

Solution: 1079j^l4°C (1975+25°F) for 4 hours; air cool 

Precipitation: 899+14°C (1659+25°F) for 10 hours; air cool. 

3.2 Specimen Fabrication 


The gage section geometries of the monotonic tensile/creep, isothermal 
cyclic, thermomechanical cyclic, and isothermal biaxial test specimens were 
chosen to correspond to the gage section geometries of the corresponding 
specimens in the concurrent life prediction contract. Major differences, 
however, in the specimen gripping systems between SwRI and PWA resulted in 
different overall specimen and casting configurations. Figure 3.2. The 
use of col letted grips at SwRI considerably simplified the specimen design, 
but added to the length and diameter of each casting. A number of casting 
trials and specimen design changes were evaluated before suitable specimen 
and casting configurations were chosen to satisfy grain size and integrity 
goals. Additional details on the fabrication are provided below for each 
specimen type: 

• Tensile and Creep Specimens - The increase in cast bar length 
from a standard configuration of 4 inches to 7.2 inches did 
not appreaciably alter specimen integrity or grain size. 

Twenty-two test specimens were prepared from the cast bars by 
centerless grinding. All specimens were subsequently electro- 
polished prior to shipment to SwRI for testing. 

• Thermomechanical Constitutive Specimen - A common, constant 
OD 2.74 cm (1.08 inch) cast bar configuration of 20 cm 

(8 inch) length was proposed from which each type of uniaxial 
constitutive specimen could be machined. Grain size and 
porosity goals could not be satisfied simultaneously with a 
casting of this configuration. Casting of "net shape" bars 
of two lengths was then attempted. Microstructural evaluation 
indicated that an adequate casting could be obtained if the 
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Modified casting procedure 12X 

FIGURE 3.1. COMPARISON OF CAST BAR POROSITY. 
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total cast bar length were limited to 6-1/2 inches. 

Figure 3.3 indicates the acceptable level of micro- 
porosity in the gage section of the 6-1/2 inch speci- 
men. Twenty- five castings of this configuration are 
available. Cylindrical buttons of .8 inch length are 
being electron beam welded to the ends of each speci- 
men to increase the stock length to the required 
8 inches. 

• Isothermal Constitutive Specimen - Difficulties in 
casting the 2.74 cm OD x 20 cm L (1.08 inch x 8.25 inch L) 
solid bars prompted a reappraisal of the recommended speci- 
men configuration. A reduction in specimen length to 
17.78 cm (7.0 inches) and diameter to 1.43 cm (.563 inches) 
was agreed upon. This change permitted use of the tensile/ 
creep specimen casting configuration. Twenty-eight cast- 
ings have been obtained and are being machined. 

• Multiaxial Tension-Torsion Specimen - A constant OD 30.2 cm 
(1.19 inch) cast bar configuration of 20.3 cm (8 inch) 
length was initially evaluated. Grain size and porosity 
levels could not be satisfied simultaneously. A second 
casting attempt was made with a revised gating scheme, but 
microporosity and grain size results were equally unsatis- 
factory. A decision was reached to cast these specimens 
with a hollow core. Hitchiner Mfg. of Milford, N.H, has 
been contracted to provide the castings. 
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(a) Gage section (b) Shoulder section 

X ; . ' % ‘f 



(c) Grip section 



FIGURE 3.3. THE LEVEL OF MICROPOROSITY AT VARIOUS POSITIONS OF A TYPICAL SPECIMEN. 
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4.0 TASK C. UNIAXIAL EVALUATION OF CONSTITUTIVE MODELS 


4.1 Tensile, Creep and Cyclic Property Testing 

The B1 900+Hf specimen requirement for the base program is presented 
in Table 4-1. A total of 39 smooth uniaxial bar specimens will be fabri- 
cated to provide constitutive response data for the loading conditions 
in Task C. A total of 20 biaxial specimens will be prepared for response 
data under multiaxial stress states for the effort under Task E. Ten 
notched benchmark specimens will be available for the work in Task F. 

All testing is being performed at SwRI. The specimens are being fabri- 
cated by PWA as described in the preceding Section 3. 

The matrices for tensile and isothermal cyclic testing are given in 
Tables 4. II and 4. Ill, respectively. All the tensile and five creep tests 
have been completed. The cyclic testing is just being initiated. Cor- 
relations with cyclic data given below are based on results from PWA 
Contract NAS3-23288. 

All the uniaxial tests are being conducted in a closed-loop, servo- 
hydraulic test machine under strain control. Figure 4.1 shows the con- 
figuration of the test specimen, the hydraulic collet grips, the induction 
heating coil and the externally mounted extensometer with quartz reach- 
rods attached to the specimen. The plastic tubing carries cooling water 
to the hydraulic collets and copper face plates. Temperature variations 
over the specimen gage section are within ASTM specifications for short 
term tests. 

4.2 Comparison of Experiment and Bodner-Partom Model 

The material constants in the Bodner-Partom elastic-viscoplastic 
constitutive theory were evaluated for B1 900+Hf. Most of the constants 
were evaluated from uniaxial tensile and isothermal cyclic data. Creep 
data were occasionally used for evaluating the constants in the thermal 
recovery terms when slow strain rate (£< 4 x 10“ 6 sec _1 ) tensile data were 
not available. A summary of the Bodner-Partom theory and an outline of 
a systematic procedure for evaluating the Bodner-Partom model constants 
analytically are shown in Figures 4.2 and 4.3, respectively, while the 
various material constants are shown as a function of temperature in 
Figures 4.4(a) and 4.4(b). It is worthy to note that for B1900+Hf most 
of the material constants in the Bodner-Partom theory are independent 
of temperature; in the temperature range of 25-1093°C, only three param- 
eters (n, Zp, and A) vary with temperature. Another material input to 
the Bodner-Partom theory is the elastic modulus which was determined 
from uniaxial tensile data and found to decrease with temperature also, 
as shown in Figure 4.5. 
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TABLE 4.1 


B1 900+Hf SPECIMEN FABRICATION REQUIREMENTS FOR THE BASE PROGRAM 


Specimen Type 

Base Program 

Task 

Fig. 3.2 

Tensile Specimens 

11 

C 

(a) 

Creep Specimens 

5 

C 

(a) 

Isothermal Cyclic Constitutive Specimens 

20 

C 

(b) 

Thermomechanical Constitutive Specimens 

3 

c 

(c) 

Multiaxial Cyclic Constitutive Specimens 

20 

E 

(d) 

Benchmark Notch Specimens 

10 

F 

★ 


* To be determined. 
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TABLE 4. II 


INTEGRATED TENSILE TEST MATRIX 


Strain Rate 



Temperature (°F) 



(sec“l ) 

RT 

1200 

1400 

1600 

1800 

2000 

4 x 10“ 7 



★ 



* 

4 x 10" 6 



★ 

1 

1 

1 

l n 
o 

X 

U* 

1 

1,* 

2 

2 

2,* 

X 

o 

1 



★ 

1 

1 

1 

4 x 10" 3 



★ 



* 

4 x 10 -2 



★ 



* 


1,2 = Test requirements, Contract NAS3-23288. 

* = Additional test conditions under this contract. 
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TABLE 4. Ill 


ISOTHERMAL CYCLIC CONSTITUTIVE TESTS 


Test 9 

Temp( a F) 


R 

^(sec’ 1 ) 

AC ] 1 d £2 

u 3 

Ac 4 

0 

effect of Temperature 

and 

Strain 

Rate 




1 


RT 


0 

4x1 O ' 4 

X X 

X 

X 

2 


800 


0 

4x 1 0 ” 4 

X X 

X 

X 

3 


300 


0 

4x1 0" 2 * 

X X 

X 

X 

4 


800 


0 

4x1 O ' 6 

X X 

X 

X 



1200 


Q 

4x1 0 * 4 

X X 

X 

X 

5 


1200 


0 

4x10"^* 

X X 

X 

X 

7 


1200 


0 

4x1 0 " 6 

X X 

X 

X 

3 


1400 


-1 i 

4x1 O' 2 * 

X X 

X 

X 

9 


1600 


- 1 

Thermal memory 
4x10-4 

erasure cycle 

X X 

X 

X 

10 


1800 


-i ! 

1 Thermal memory 

erasure cycle 



n 


2000 


-i 

4x1 0 * 6 

X X 

X 

X 

0 

Effect 

of Mean Strain 







12 


1200 


-i 

4x1 0 * 4 

Same as 5 except 

fully reversed 

13 


1400 


0 

4x10-2* 

Same as 8 except 

one way 

tension 

14 


1400 


— 

4x10-2* 

Same as 8 except 

one way 

compression 

15 


1600 


— 

4x1 0 “ 4 

Same as 9 except one way 

ccmcression 

15 


2000 


«<a 

4x10-6 

Same as 11 except 

one way ccmorsssion 

0 

Strain 

Holds to Probe 

for 

Sack Stress 




17 


1200 


0 

4 X 1 Q— 4 

X X 

X 

X 

13 


1400 



4x1 0 “ 4 

X X 

A 

X 

13 


1600 


-1 

4xl0 " 4 

X X 

’/ 

A 

V 

A 

20 


2CCQ 


-1 

ixlO "* 1 

X X 

X 

V 

A 


*It may be desirable to recuce this uoDer strain rate (4x10”” S" 1 ) based on 
actual component rates and on measurement accuracy for small inelastic strains. 
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FIGURE 4.1. TENSILE SPECIMEN UNDER TEST AT 1600°F. 
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BODNER - PARTOM'S MODEL 


Flow Law: 


e.. = £. 6 + 

1J 1J 


"1J 


6 1j = X S ij ; *kk = 0 

With S ij = °ij ' I 5 ij°kk 
Kinetic Equation: 



D p = Ii p ; p 
U 2 2 e ij e ij 


J 9 =7$, .S. . 

2 2 lj ij 

2 P 

\ = d 2 /j 2 


Evolution Equations of Internal Variables: 
a. Isotropic Hardening 

,1 


.1 

Z = m 


V* Z 3 - Z ‘]»p - 


r-rl 


A 1 Z 1 


V - 1 2 

T, 


wh ere d = (1 - a) W p sin0 


9 

= cos"' 

(v.,v 

ij> 

or 9 

= cos"' 

(u. • 

u. .) 



ij 





v ij 

= B id /(e 

kl B kl 

> 1/2 , 

* 7 ij 

= v< 8 

kl 8 kl 

) 1 /2 

u ij 

= o 1d /(o 

kl a kl 

) 1/2 , 

1 U ij 

= ^/(a 

kl d kl 

} l/2 

Z 1 1 

O 

II 

rvi 

0 

; *P 

= °1j 

P 

*U ; 

w p ( 0 ) 

= 0; 

a(0) 


b. Directional Hardening 


®ij = m 2 (Z 3 u ij - 8 i-j> “p - A 2 Z 1 


(BklBkl) 


1/2 


ij 


with Z D = 8 i jU . j ; Z D (o) = 0 , 6^(0) = 0 

m ? r ^ 

m 2 = “ 2 “ (1 + sin0) |J + expt-m^Z^J 
Material Constants: D Q> Z Q , Z] , Z 2 * Z3, m 1 , m 2 , m3 


A], A^, r-j , r 2 , n, and elastic constants 


FIGURE 4.2. A SUMMARY OF BODNER-PARTOM' S MODEL. 
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FIGURE 4.3. A FLOW CHART SHOWING A SYSTEMATIC PROCEDURE FOR 
EVALUATING THE MATERIAL CONSTANTS IN THE 
BODNER-PARTOM THEORY. 
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MPa 



FIGURE 4.4(a). THE VALUES OF Z 0 , Zi , Z 3 , mi, AND m 2 
FOR B1 900+Hf AT VARIOUS TEMPERATURES. 



FIGURE 4.4(b). THE VALUES OF D 0> n. A, AND r FOR 
B1 900+Hf AT VARIOUS TEMPERATURES. 
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The Bodner-Partom model was used to simulate the uniaxial tensile 
behavior of B1 900+Hf as a function of both temperature and strain rate. 
Stress-strain curves were first calculated and the results were then used 
to obtain the yield stress" at .2% plastic strain for various imposed 
strain rates and temperatures. Figure 4.5 shows that the calculated 
.2% eP yield stresses at e = 8.3 x 10"' 3 sec" 1 at various temperatures are 
in good agreement with the experimental values. The effects of strain 
rate on the .2% cP yield stress at 760°C, 871°C, and 982°C are summarized 
in Figure 4.6. The results indicate that the rate effect is relatively 
small compared to that of temperature, in good agreement with experimental 
observations. Figure 4.7 shows the stress-strain curves at e = 8.3 x 10"^ 
sec" at 538, 871, and 982°C. Both the experimental and the calculated 
curves show that the flow stress at a given strain decreases with tem- 
perature while the apparent hardening rate is increased as the temperature 
is increased. 

The Bodner-Partom model was used to calculate the isothermal cyclic 
hysteresis loops at various strain ranges using the same set of constants 
shown in Figures 4.4(a) and 4.4(b). The calculated cycl ical ly-saturated 
hysteresis loops at various strain ranges are compared with experimental 
data in Figures 4.8 and 4.9 for 538°C and 871°C, respectively. The com- 
parisons indicate reasonably good agreement between theory and experiment 
for both small (near yield) and large plastic strain hysteresis loops. 

Cyclic hardening is also predicted by the Bodner-Partom theory. Figure 4.10(a) 
shows that both the calculated hysteresis loops after 1-2 cycles and at cyclic 
saturation; these calculated loops are in fair agreement with the experimental 
loops shown in Figure 4.10(b). 

4.3 Comparison of Experiment with Walker Model 

A very large quantity of constitutive information has already been 
generated as part of the data acquisition efforts in the HOST Isotropic 
Fatigue Contract (NAS3-23288) . Much more data will be generated as that con- 
tract continues, as well as in the present contract. An important need in 
the present contract effort and in other related HOST programs is the ability 
to deal expeditiously with this information while still maintaining a high 
level of insight into the underlying physical processes. A computer program 
to do this utilizes the C0PES/C0NMIN Optimization Code written by L. E. 

Masden and G. N. Vanderplaats. The C0PES/C0NMIN deck has been coupled with 
PWA codes which provide interactive access for real time use of the COPES/ 
CONMIN system and direct access to PWA graphics/file manipulation routines 
and the data base established for the HOST Isotropic Fatigue Contract data. 

A flow chart describing these decks and their interrelations is shown in 
Figure 4.11. A summary of Walker's model is shown in Figure 4.12. 

In addition to the motivation to use optimization techniques resulting 
from the large volume of constitutive data available and the heuristic nature 
of the current constitutive theories, the limited deformation observed in the 
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4.5 THE ELASTIC MODULUS AND THE .2% e P YIELD STRESS AS A 
FUNCTION OF TEMPERATURE. The calculated . 2 % eP 
yield stresses using Bodner-Partom theory are in 
good agreement with experimental data. 



OF STRAIN RATE FOR B1900+Hf AT 760, 871 and 982°C. 
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Stress, MPa 



FIGURE 4.7 A COMPARISON OF THE CALCULATED (BODNER-PARTOM THEORY) 
AND THE EXPERIMENTAL STRESS-STRAIN CURVES OF 
B1 900+Hf AT 538, 871 AND 982°C. 
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STRESS. MPR Stress, MPa 



(a) Experiment 



Strain, % 

(b) Bodner-Partom Theory 


FIGURE 4.8. A COMPARISON OF THE CALCULATED (BODNER-PARTOM THEORY) 
AND THE EXPERIMENTAL SATURATED HYSTERESIS LOOPS FOR 
B1900+Hf AT 538°C. 
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(a) Experiment 


Bodner-Partom Theory 
B1900+Hf 


i = +8,3 x 10 4 sec~^ 


+ .4% 

+ .35% - 



Strain, % 


r •* i 


(b) Bodner-Partom Theory 


FIGURE 4.9. A COMPARISON OF THE CALCULATED (BODNER-PARTOM THEORY) 
AND THE EXPERIMENTAL SATURATED HYSTERESIS LOOPS FOR 
B1 900+Hf AT 871 °C. 





Stress, MPa 



(a) Experiment 



(b) Bodner-Partom Theory 


FIGURE 4.10. A COMPARISON OF THE CALCULATED (BODNER-PARTOM THEORY) 
AND THE EXPERIMENTAL HYSTERESIS LOOPS AFTER 1-2 CYCLES 
AND AT CYCLIC SATURATION. 
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FIGURE 4.11(a). A FLOW CHART SHOWING THE COPES/CONMIN 
OPTIMIZATION CODE AND ITS INTERACTIONS 
WITH PWA GRAPHICS/FILE MANIPULATION 
ROUTINES AND THE DATA BASE. 



FIGURE 4.11(b). A FLOW CHART SHOWING THE OPTIMIZATION PROCEDURES 
UTILIZED IN SELECTING MATERIAL CONSTANTS IN 
WALKER'S MODEL. 
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oc: 


WALKER'S MODEL 



Material Constants : 


o 

A > p, ^ i n , m, n -J } ^ 2 * n 3 » » ^5* 

K-j , K^, depend on temperature. 


n 7’ n 8’ 


FIGURE 4.12. A SUMMARY OF WALKER'S MODEL. 
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testing of the precipitation aged materials is also a prime factor No 

sire's S etT) ?Li e »ni C ?3 b ^.f r f or ™ ed (^ess <»P tests to probe'for back 
_x ’ . e + c ’' that will identify in a straightforward manner the evolution 

Of the internal variables for materials of this type. If "trial fnd error" 

is required, optimization techniques can considerably facilitate the search. 

11C . As a of the studies performed to date, several aspects of the 

severll^aHati^rnf^h t0 be especially important. In addition, 

several variations of the basic approach were evaluated and still other 

variations identified for continued investigation. A description of these 
findings and the procedures used follows. P ,on 0T tnese 

The interactive computer code which was developed as part of this 
effort is outlined in Figure 4.11(a) with additional details presented in 
Figure 4.11(b). The procedure normally followed in the use of the code begins 
with the selection of the data set from the HOST data base to be used in 
the optimization process. This data set is returned to the interactive 
control program as shown in Figure 4.11(a). The user then exercises a series 
of options as outlined in Figure 4.11(b). These include: 1) the option to 
reduce the number of data points in the file selected in order to control 
computation time, 2) the option to provide a weighting matrix to bias the 
least squares fit toward those aspects of the data (data fit at strain 

limftc 3 p01nt ? for exam P 1e ) deemed most significant, 3) the ability to set 
formation ^ ons J ra ] nt ^ ° n ^ he magnitude of the constants based on prior in- 
[nforma??on ,nrh S1C and finally ’ 4) to sup P ly 9eneral control 

foSion desired C ° nVergence cnterion a " d the type of optimization in- 

„ nMrr The deri va tion of constitutive theories stems generally from two 
sources. The first source is the isolated mechanical behaviors observed 
in laboratory tests. Second, constitutive theories are derived from con- 

enti rp^ nd h ?° n6S °l P bysical metallurgy. Neither of these sources is 
entirely exclusive of the other or clearly sufficient in itself. 

„ j i It: I s f or ^ h , ls reason that a viable approach to the determination of 
f constants might proceed on the notion that the most appropriate tests 
from which to determine constants are those which most closely represent the 
application where the model is to be utilized. Preferably, this set of base 

condit?nn^ C fh t f StS f ° r mode i constant determination would span the range of 
conditions that are reasonably expected in hot section blade and vane ap- 

equipment S ’ 1mited 0n y by the usua1 1 imi ts imposed by use of laboratory 

Based on the considerations outlined above, a study was initiated to 
determine to what extent useful constitutive model constants could be de- 

thlt^ho mn 1 ? 9 C ^ C ^ data al0ne> More specifically, it was hypothesized 
t t re | able t 0u ^ c ® °f constants would be obtained from a fatigue 

test with the most involved loading wave form. This approach was tested 

rpJnt 6 °°h F and 1 18 ?°° F data from the HOST Isotropic Fatigue program. The 
results and conclusions are presented below. H y 
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The response of the material during a fully reversed dwell fatigue 
test (Specimen 27B) was selected as the base or reference data for selection 
of the model constants at 1600°F. 

The model constants were found by optimizing a sum of squares param- 
eter based on the residual deviations between the model predicted stress 
and the measured stress at each recorded data point. No constraints on the 
Walker constitutive model constants were imposed. The optimum fit of the 
data is compared against the measured data in Figure 4.13. The fit shown 
appears reasonable with the maximum difference at any strain between model 
predicted stress and data no more than 5% of the total stress range. However, 
a more interesting test of the validity of determining constants from cyclic 
data is to use the model constants to predict cyclic data at conditions which 
were not used to determine constants. This was done for a variety of test 
conditions, shown in Table 4. IV. The resulting comparisons of model pre- 
diction versus data at 1600°F are shown in Figures 4.14, 4.15, and 4.16 for 
one way rapid cycle, one way slow cycle, and fully reversed rapid cycle large 
strain conditions. The maximum deviations between model predicted stress 
and data are similar in magnitude to that found for the reference Specimen 
(27B) used to determine model constants at 1600°F. 

While the approach of determining constants from cyclic data was suc- 
cessful in predicting the behavior at other cyclic load conditions at the 
same temperature, it was not successful in predicting the material behavior 
during a mono tonic tensile test. Comparisons of the ability to predict high 
strain tensile data at two strain ramp rates from cyclically determined con- 
stants are shown in Figures 4.17 and 4.18. As readily noted, the high strain 
data fit is poor. Figure 4.19 presents a repeat of the condition of Figure 
4.17, but with model constants determined from the data itself to illustrate 
that the optimization can readily tailor constants to individual data sets. 

An additional approach to evaluating the validity of a constitutive 
model formulation is to compare the optimized constants from each test con- 
dition and observe the scatter between individual loading conditions. This 
comparison is shown for a partial list of the constants from the cyclic con- 
ditions at 1 600°F in Table 4.V. It was found that the optimum constants 
were not influenced by varying cyclic strain rate, R-ratio, and loading 
wave form (triangle wave vs. hold periods). 

Not unexpectedly the strain hardening exponent and the principal com- 
ponent of the drag stress, K-j , were influenced by a change in strain range 
as shown by a comparison of the constants for Specimens 31 B and 27B in 
Table 4.V. This trend is similar to that observed from the comparisons of 
model prediction versus data for the tensile tests. 

Constants for Walker's model were also determined at 1800°F. A com- 
parison of the model using optimized constants with the corresponding data 
is shown in Figure 4.20 for a slow cycle fully reversed test. Figure 4.21 
shows a comparison of model prediction with the tensile test data at 1800°F, 
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TABLE 4. IV 

SPECIMEN TEST CONDITIONS USED IN WALKER'S MODEL EVALUATION 
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Cycle 

Tensile, Rapid 41 B - - -005 in/in/min 



TABLE 4. V 


PARTIAL LIST OF CONSTANTS FOR WALKER'S MODEL 
DETERMINED BY OPTIMIZING ON CYCLIC DATA FROM INDIVIDUAL TESTS 


SPECIMEN TEST DESIGNATION* 




27B ( 1600 F ) 

21B (1600°n 

T1A (1600 8 n 

~m (1600 a F) 

29B ( 1800°F) 

n 

(Strain Hardening 
Exponent) 

6.203 

6.250 

6.090 

7.107 

4.476 


Initial Equilibrium 
Stress Value (psi) 

-2828.7 

-2829.1 

-2828.7 

-2827.0 

-3435. 

K i 

Principal Component 
of Drag Stress 

.2288 x 10 6 

.2286 x 10 6 

.2288 x 10 6 

.1704 x 10 6 

.1839 x 10 6 

n i 


.01340 

.01337 

.01338 

.01336 

.01338 

n 2 


.1018 x 10 9 

.1019 x 10 9 

.1018 x 10 9 

.1015 x 10 9 

.9890 x 10 8 

n 3 


.1990 x 10 4 

.1990 x 10 4 

.1990 x 10 4 

.1997 x 10 4 

.2022 x 10 4 

n 6 


.472 

.472 

.472 

.472 

.472 

n 7 


10.215 

10.219 

10.215 

10.202 

10.154 

n 9 


340.65 

340.26 

340.65 

240.0 

333.4 

n ll 


.9989 x 10 7 

.1000 x 10 8 

.9989 x 10 7 

.9975 x 10 7 

.9973 x 10 7 

m 


1.0224 

1.0197 

1.0220 

1.0337 

1.2000 

E 

Young's Modulus 

20.13 x 10 6 

20.04 x 10 6 

20.12 x 10 6 

20.13 x 10 6 

16.84 x 10 6 


Remaining constants were not included either because they have been found not to 
contribute for precipitation hardened alloys or pertain only for non-proportional 
loading. K K 

*Test conditions for each specimen shown in Table 4. III. 

Constants n . , n n~, n., n^, n,, and m influence the evolution of the equilibrium 
stress in Walker's modeT n 7 controls the rate of cyclic hardening. 
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FIGURE 4.13. OPTIMUM FIT OF 1600°F FULLY REVERSED COMPRESSION 
DWELL DATA USING WALKER'S MODEL (SPECIMEN 27B). 



FIGURE 4.14. FIT OF 1600°F ONE WAY TENSILE, RAPID CYCLE 
DATA (SPECIMEN 21B) USING WALKER'S MODEL 
WITH CONSTANTS DERIVED FROM SPECIMEN 27B. 
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TIME, SECONDS 


FIGURE 4.15. FIT OF 1600°F, ONE WAY TENSILE, SLOW CYCLE DATA 
(SPECIMEN 21 A) USING WALKER'S MODEL WITH 
CONSTANTS DERIVED FROM SPECIMEN 27B. 


S 


V, 


TIME, SECONDS 

FIGURE 4.16. FIT OF 1600°F, FULLY REVERSED, RAPID CYCLE DATA 
AT +.4% (SPECIMEN 31 B) USING WALKER'S MODEL WITH 
CONSTANTS DERIVED FROM SPECIMEN 27B. 
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FIGURE 4.17. FIT OF 1600°F MONOTONIC TENSILE DATA AT 
A STRAIN RATE OF .5% PER MINUTE USING 
WALKER'S MODEL WITH CONSTANTS DERIVED 
FROM SPECIMEN 27B. 



STRAIN, IN/IN 

FIGURE 4.18. FIT OF 1600°F MONOTONIC TENSILE DATA AT 
A STRAIN RATE OF .05% PER MINUTE USING 
WALKER'S MODEL WITH CONSTANTS DERIVED 
FROM SPECIMEN 27B. 
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DATA, x 10 4 PSI 


1.0 



STRAIN, IN/IN 

FIGURE 4.19. OPTIMUM FIT OF 1600°F MONOTONIC TENSILE DATA 
AT A STRAIN RATE OF .5% PER MINUTE USING 
WALKER'S MODEL. 



FIGURE 4.20. OPTIMUM FIT OF 1800°F FULLY REVERSED, SLOW CYCLE 
DATA USING WALKER'S MODEL (SPECIMEN 29B). 
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FIGURE 4.21. FIT OF 1800°F MONOTONIC TENSILE DATA AT A STRAIN RATE 
OF . 5 % PER MINUTE USING WALKER'S MODEL WITH CONSTANTS 
DERIVED FROM SPECIMEN 29B. 


the constants derived from the cyclic test were used for this comparison. 

The same trends observed at 1600°F were found at 1800°F. 

The method just described utilizes optimization of an objective 
function based upon the least squares parameter calculated from the residual 
deviation at each measured data point. Four rather specific strategies for 
the use of numerical optimization were tried. 1) Unrestrained optimizations 
on all constants simultaneously. 2) Progressive optimization on several 
types of cyclic tests using increasingly tight constraints on those constants 
determinable from prior data fitting, e.g., elastic constants. 3) A pro- 
cedure which allowed constants to be fixed (equality type constraints), once 
determined, so that optimization can proceed on the remainder of the con- 
stant set. Of these three approaches, the first has surprisingly proven to 
be the most advantageous to date. The second approach was initially expected 
to be the ideal one, but the optimization algorithms used in this study tended 
to terminate the search procedures prematurely when the search limits (upper 
and lower bounds) were set close to one another. This limitation of the 
optimization procedures has been noted by the author of the CONMIN code, Dr. 

G. H. Vanderplaats . 

The third strategy involved the temporary removal of the constrained 
variable from the list of those to be optimized. Removing a model constant 
from the variable list this way amounts to a practical way of providing an 
equality constraint on the constant. This approach yielded good results 
initially and will be explored further. The principal disadvantage of this 
approach is that it assumes that a particular constant is known exactly rather 
than within a narrow uncertainty band. This should not be a serious drawback 
in practice, however. 

A fourth approach was the use of weighting factors applied to those 
points in the data cycle deemed most significant. The trials (e.g., increased 
significance given to data fit near strain reversal or null load points) using 
this approach showed no advantage on the limited data set evaluated. 

Other strategies which remain to be explored are the use of constraints 
on dependent variables associated with the model and model constants and 
secondly, the use of alternative forms of the objective function. Two ex- 
amples of dependent variable constraints which may prove useful are: 1) a 

narrow band around the actual measured response (e.g., measured stress) which 
might be made progressively more narrow in order to guide the optimization in 
the desired direction, and 2) a similar bound on the maximum residual deviation 
(e.g., the maximum absolute value of the differences between the measured and 
calculated stresses for all data points). This same maximum deviation should 
also be explored as an alternative to the sum of the squared deviations as an 
objective function to be minimized. It is considered likely that this would 
be a more sensitive parameter than the more usual sum of squared residuals. 

It may have the drawback, however, that it is overly sensitive to random 
fluctuations, or noise, in the measured stress data set and may not converge 
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due to minor discontinuities. A second alternative to the objective function 
which could possibly avoid sensitivity to data set noise would be weighted 
combination of the sum of squared residuals and the maximum deviation. 

To summarize the experience with the use of optimization for deter- 
mination of constitutive model constants, it is to a large degree a very 
successful approach. It is too early to say that the specific strategies 
used are the best ones to use. No one method has been found completely 
satisfactory at this time. However, it is clear that reasonable constant 
values can be determined with the relatively small expenditure of analysis 
time, which was the primary goal of the study to date. 
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5.0 TASK D. 


IMPLEMENTATION IN FINITE ELEMENT COMPUTER CODE 


. , , An ^ constitutive model being developed must be included in a struc- 

IZ 1abp 3 yS ! S prog ^ m : In this and tw o previous contract efforts [1 and 21 

the MARC nonlinear finite element computer program was the vehicle for inrnr’ 

ZleZlT- V ]) C a° P JrS C d m0del " Th ? ™ lowi *" g -ct?on Siff VtilZe'Z'o 

cnree parts. 1) a bnef discussion of the MARC code and the procedure fnr 
implementing a viscoplastic model, 2) a review of various integrating methods 

rnHp V1 ?i° P astlc the ° nes considered to date in conjunction with the MARC 
code, 3) a summary of the approach being taken in Task D. 

5 - ] Description of the MARC Program 

In References 1 and 2, the viscoplastic constitutive theories were in 
corpora ted into MARC program by means of an initial stress technique All of 

an in? tfaT LT IT**. 1 ! the COnstitutive equations is incorporated into 
an initial load vector and treated as a pseudo body force in the finite elp- 

for a 6 “st the v1 > c0 P ,a5tfc constitutive theories 

incrMentirennstitut?I» d l ff ^" tla equat,ons - U is necessary to form the 
incremental constitutive equation appropriate to the finite element load 

increment by means of a subincrement technique. 

cn ,.. the subincrement technique the finite element load increment is 
split into a number of equal subincrements and the viscoplastic constitutive 

presentation TZT* ° Ver Sm "■ *"M"«=r-e«ts to ?om an tccllllt ll- 

?^H e " tatl ° t ? e 1ncremental constitutive equation over the finite element 
nZtl nC ?T n l- Integration °ver each subincrement is accomplished by a 
that 6 th 0f -? Ue i' Provided th e subincrements are sufficiently small (so 
tSrL^ he s ^ abl J lt ^ ]. eve l of the integration method is not exceeded), the 
technique has been found to work efficiently and accurately, even for large 

nird^fr- 6016 "! °K d lncrements * However, it is difficult for the user to 9 

use as f fIS 1 c l1 M SUblnCr ? nenta1 steps ’ and there is a considerable incentive to 
! f subincrements as possible, consistent with the stability of the 

differential equations comprising the constitutive theory. 

relationshins R intn d thf ll0WS the . user to im Plement very general constitutive 
in thic h PS l" t0 pro 9 ram b y means of the user subroutine HYPELA. With- 
D • aPd ! n %- he ? S6r mUSt specif y values of the elasticity matrix 

const? tutive 1 re1ationship 6SS ’" C, ' eraent V6Ct ° r 1n the incremental ™ctor 


Aor. = D.^A . - dj-aAe) - A^ 


(5.1) 


The inelastic stress increment vector Ac- is computed in HYPELA from the visco- 
plastic constitutive relationships. 1 


In Equation (5.) a denotes the coefficient of thermal 
vector Kronecker delta symbol, 


expansion and <$• is the 
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^ 1 if Os j s 3 
0 if 3 < j s 6 


(5.2) 


For the class of nonlinear viscoplastic constitutive relationships 
under consideration in this contract, the incremental inelastic stress 
vector ACj depends in a highly nonlinear manner on the incremental strain 

vector Ae-j . 

The solution of the incremental equilibrium is accomplished within 
the MARC code by the following algorithm. At the start of the increment 
the user subroutine HYPELA is entered to determine the elasticity matrix 
D-: i and the incremental inelastic stress vector Aq. On entry to the sub- 
routine the input consists of the strain increment vector A ei , the 
ture increment A©, the time increment At over which the incremental e *^' 
nal load vector APi is applied to the structure, and the values of the stress, 
strain, temperature and viscoplastic state variables at the beginning of the 
increment. Since the incremental strain vector, Aei = B ii 
accurately determined after the solution to the incremental equilibrium - 
lationship has yielded the correct incremental solution AUj , the strain in- 
crement vector A £i initially used to generate the A^ is assumed to be the 
value obtained for Ae,* in the preceding increment. On exit from subroutine 
HYPELA the elasticity matrix D-j j and the estimated inelastic stress ' , ncre- 
ment vector Ac- are passed into the main program. After the values otUij 
and Aci are obtained for each integration point in the structure, the incre- 
mental J equil ibrium relationship is assembled and solved for the incremental 
node displacement vector AUi . The incremental strain vector, Aei = B ij A “j> 
is then computed and compared with the initial guess for Aei used to ^nera 
the inelastic incremental stress vector A^,-. If this incremental strain vec 
tor is equal, within a user specified tolerance, to the incremental strai 
vector used to compute Acj in the assembly phase, the solution 1S assumed to 
have converged. Otherwise, the updated strain increment vector, .obtained from 
the solution of the equilibrium relations is passed into subroutine HYPELA, 
new vector, acj , is computed, and the equilibrium equations resolved to yield 
an improved value of auj and Aei . The Process is repeated ° f 

the vector Aei on the assembly phase is equal, within a user specified toler- 
ance, to the value of the vector Ae n - on the solution phase. After c °™ergenc 
is achieved, the temperature, stress vector, strain vector and vlsc ?P^ a stic 
state variables are updated by adding the incremental values generated during 
the current increment to the values of these variables at the beginning 
increment. The program then passes on to the next load increment where the 

process is repeated. 

5.2 Integration Methods for Viscoplastic Theories in the MARC Co(te 

The values of D n * j and Aq in the incremental constitutive relation. 

Act. = D . . ( Ae . - 6. aAG) - A£. 

1 1J J J 1 
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are obtained by a subincrement method. Incremental values of the variables 
a 0 , At and (Ae-j and 6-j aAO) for the current finite element load increment 
are split into NSPLIT equal values, and the constitutive equations are inte- 
grated over the NSPLIT subincrements to provide accurate values of D-j j and 
A^j. Each load increment in a MARC analysis is divided into NSPLIT subincre- 
ments. The integration of the constitutive equations is currently performed 
by using explicit Euler forward differences with a step size determined by 
dividing the MARC Load increment by NSPLIT. 

Three methods for integrating the constitutive equations were examined 
in [2]: (1) an explicit Euler difference scheme with error estimates for re- 

vising the time step, (2) a backward difference scheme, and (3) integration of 
an integral form of the equations. 

The forward difference integration, similar to [3], is based on an 
error estimate, e, given by 


^ 3AJ 2 
= AR + \ 


(5.4) 


where 


AR = J Ae.P Ae.P 

* 3 ij ij 


p 

is the inelastic strain 
S^j is the deviatoric stress 

A of a quantity is the change in the quantity over time step At. If the 
error estimate is too large 


e > e ] (5.5) 

then the time step is halved and the step is repeated. If the error estimate 
is too small 


e < 62 


(5.6) 


then the time step is doubled for the next integration step. 
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A backward difference algorithm was also developed in [21. Although 
the forward difference worked adequately for all test cases considered, the 
backward difference convergence was slow in cases where the term ng, in 
Walker's equation, was not equal to zero and the strain rate was smaller 
(i.e., on the order of 10' 6 /sec). A quadratic Newton's method was also used 
to solve for aG and but there was no benefit over a linear Newton's 

method. Therefore, in these cases, a sufficiently small step size was re- 
quired. 


The equations for the modified Walker's theory have also been recast 
in integral form, [1]. 

The implicit integration for the evaluation of the integral form of 
the equations can be illustrated by considering the first order differential 
equation 


y(t) = f Q (t) - y ( t ) f 1 ( t ) (5.7) 

where f Q and f 1 are known functions of time, subject to the initial condition 


y(o) = y 0 


(5.8) 


Equations (5.7) and (5.8) have the solution 


y(t) - y 0 e-«<> ♦ / .-tQUMXO] (t)d{ 

0 

where 

Q(n = f V T)d 

o 

From [11, for large times. Equation (5.9) is 


y(t) « -2 
f l 




f l f 0 + f l 


(5.9) 


(5.10) 


(5.11) 
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Equation (5.9) can be readily cast in an incremental form as 


^+1 


-CQ(t, lAl )-Q(i)j 


~ t Q( , 1 } “Q( ttj ) ] f “ L V\ Uij,-. 

y (t N+1 ) -y(t N )e *j e f o ,()d< (5-12) 


U N 


An approximate numerical integration of Equation (5.12) using Equation 
(5.11) is 


■^N+l ^ 


[Q(t N+ i )-Q(t w )] fnftu^) 


0 1 H+l 

VW 


(5.13) 


where 


- e 


- w<W-Q<V ] W , . Af o . Af i 

f 1 ( t N } f 0 (t N^ W 


Af 0 " ^O^N+l^ " and = 


.The evaluation of Q(t) is given by Equation (5.10) 


Q(t N+l ] 


Q(t N ) + 


rh+i 

l f ’ 




(5.14) 


and Equation (5.12) can be approximated by using a convenient numerical inte- 
gration rule. For example, a rectangular rule gives 


^Sj+1 ^ + (t n ) At 


(5.15) 


where 


At = 



t 


hi 


or if desired, a trapezoidal rule can be used, then 
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(5.16) 


Q(t N+l ) = Q( V + 7 C f !<t N ) + ^ (t N+1 ) J At 
For large time steps, from Equation (5.13) 


y( W - 


f 0 (t M+l 

*1 ( ^+1 


) 

T 


(5.17) 


which is the first term in Equation (5.11), and for small time step 


y(t N+l ) ~ y( WV At " VV At ( 5>18 ) 


which is a first order Euler integration. 


5.3 Summary of Approach to be Taken 


Several numerical methods for implementing Bodner's and Walker's 
theory in the MARC finite element computer program are being considered in 
this Task. These include, but may not be limited to: 

1) Explicit Euler integration with both a fixed and self- 
adaptive time step. This approach has been used success- 
fully in the past with Walker's theory. 

2) The NONSS method of Miller and Tanaka. This approach has 
been used successfully with Bodner-Partom's model. Recent 
work by Brockman have suggested a method for improving the 
efficiency of this approach [4]. 


3) Implicit integration of the integral form of the equations, 
as described above. 


4) For completeness, a second order Runge-Kutta method with a 
self-adaptive time step will be considered. The existing 
HYPELA model for Walker's theory will be modified first. 
Each theory will be finally coded with at least two numeri- 
cal integration algorithms. 
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6.0 TASK E. MULTIAXIAL EVALUATION OF CONSTITUTIVE MODELS 


The literature review conducted in Task A reveals that several struc- 
tural alloys including Hastelloy-X exhibit considerably more cyclic hardening 
when deformed under nonproportional loading paths than under proportional 
loads. To take into account this type of loading-path-proportionality- 
dependent hardening behavior, both the Bodner-Partom and the Walker models 
were modified by introducing new parameters in the evolution equation of the 
internal variables such that additional cyclic hardening would result under 
nonproportional paths. Summaries of the updated Bodner-Partom and Walker 
models are shown in Figure 4.2 and Figure 4.12, respectively. Preliminary 
calculations of cyclic stress-strain response of Hastelloy-X subjected to 
90° out-of-phase combined tension/torsion loading indicate good correlation 
between experiment and theories, as shown in Figure 6.1. However, it should 
be noted that the cyclic stress-strain curve shown in Figure 6.1(b) for the 
Bodner-Partom model was obtained using uniaxial tensile data as input only, 
while the result shown in Figure 6.1(c) was obtained by fitting the Walker 
model to the experimental data in Figure 6.1(a). 

The phase angle between the deviatoric stress and the plastic strain 
rate vectors of Hastelloy-X under 90° out-of-phase tension/torsion loading 
was also determined and the results are shown in Figure 6.2. Comparison of 
experimental data with preliminary calculations of Walker and Bodner-Partom 
models, also shown in Figure 6.2, indicate that the Walker model predicts non- 
coaxiality between the deviatoric stress and the plastic strain rate vectors, 
and the predicted phase angles which are approximately 5-15° are smaller than 
those observed experimentally. On the other hand, in the Bodner-Partom model, 
the plastic strain rate vector is always taken to be coaxial with the devia- 
toric stress vector. As a result, the predicted phase angles are always zero, 
as shown in Figure 6.2. 
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STRAIN RESPONSE OF HASTELLOY-X TESTED UNDER 90° OUT-OF-PHASE TENSION/TORSION LOADING AT 
ROOM TEMPERATURE: a) experimental data; (b) Bodner-Partom theory; (c) Walker theory. 
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FIGURE 6.2. 


A COMPARISON OF BODNER-PARTOM AND WALKER 
THEORIES WITH EXPERIMENTAL RESULT OF THE 
PHASE ANGLE BETWEEN THE DEVIATORIC STRESS 
AND THE PLASTIC STRAIN RATE VECTORS OB- 
SERVED IN HASTELLOY-X TESTED UNDER 90° 
OUT-OF-PHASE TENSION/TORSION LOADING AT 
ROOM TEMPERATURE. 
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7.0 SUMMARY OF CURRENT RESULTS 

During the first year of effort, the following tasks were accom- 
plished and preliminary results obtained: 

• A literature review of existing unified constitutive 
models for time and temperature dependent inelastic 
behavior of high temperature alloys was completed. 

All models contain three essential components - a 
flow law, a kinetic relationship between stress, 
inelastic strain rate and temperature, and evolu- 
tionary equations for the internal variable(s) de- 
scribing work hardening. Ten different models were 
compared relative to these components. 

• A review of numerical integration procedures applicable 
to these "stiff" equations was completed also. 

• The particular models of Bodner and Partom and of 
Walker were selected for further detailed study and 
implementation in a FE code. 

• The alloy selected for testing under the basic program 
was cast B1900+Hf. Specimens are being fabricated at 
PWA from a single heat of this material. 

• Uniaxial tensile and creep tests have been completed 
at SwRI for a matrix including temperatures from room 
temperature to 2000°F and strain rates from 10” 7 to 
10“ 3 sec~'. 

• Initial correlations of the uniaxial data, both 
monotonic and cyclic, with the Bodner-Partom and 
Walker models indicate reasonable to good correlation. 
Procedures for determining best-fit constants for 
these models are being developed but need further re- 
fi nement. 

• Both models are being implemented in the MARC finite 
element computer code. At least two numerical inte- 
gration algorithms will be used with each model. 

• Preliminary data from the literature indicates cyclic 
biaxial hardening to be greater than under uniaxial or 
proportional cycling. Existing models do not include 
this effect. Several modifications were attempted to 
correct this situation with reasonable success. 
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8.0 


FUTURE WORK 


During the second year of the base program, the following work will 
be completed: 

• Experimental test data will be obtained for B1900+Hf under 
isothermal and non-isothermal uniaxial strain controlled 
cyclic conditions, and proportional and non-proportional 
biaxial tension-torsion loading. 

• The Bodner-Partom and Walker models will be implemented in 
the MARC finite element code. This code will then be util- 
ized to correlate all the uniaxial and biaxial data. 

• A notched specimen geometry will be selected and tested to 
provide benchmark data under a non-homogenous strain condi- 
tion. This benchmark data will be compared with a FEM 
analysis of specimen utilizing the optimum constitutive 
model . 

• A hot section component, probably a turbine blade, will be 
analyzed using the FEM code and optimum constitutive model 
to demonstrate advanced modeling capability. 

• Finally, the FEM code with advanced unified constitutive 
model will be installed on the NASA Lewis computer facility. 
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1.0 INTRODUCTION 


Constitutive theories based on the classical concepts of plasticity and 
creep generally decompose the inelastic strain rate into a time-independent 
plastic strain rate and a time-dependent creep rate with independent constitu- 
tive relations describing plastic and creep behavior. While this approach 
can be rationalized on historical grounds and perhaps on computational conve- 
nience, experimental evidence collected on structural alloys at elevated tem- 
perature indicates inherent time-dependency and creep/plasticity interactions 
[1,2], This suggests that inelastic deformation might be primarily controlled 
by a single overall mechanism and should be treated in a unified manner. 

Although much of the essential physics of time-dependent inelastic de- 
formation of metals has been known for some time, attempts at a unified ana- 
lytical formulation have been relatively recent. This seems to have been due 
to the success of elastic-plastic analyses in many areas of engineering and 
the tendency to treat time-dependent effects as special phenomena. Probably 
the first attempt at a unified treatment was that of J. J. Gilman and W. G. 
Johnston in the late 50 # s and early 60' s who did some of the basic work in the 
field of "dislocation dynamics" [3-5]. They showed that a reasonably real- 
istic stress-strain curve could be obtained by integrating expressions for the 
elastic and plastic strain rates which were both considered to be nonzero 
throughout the loading history. 

In recent years, a number of formulations of elastic-viscoplastic con- 
stitutive equations have been presented in the engineering literature. Such 
equations are sometimes referred to as "unified" since inelastic deforma- 
tions are represented and treated by a single kinetic equation and a discrete 
set of internal variables. In this context, creep, stress relaxation, and 
plastic flow are different manifestations of time-dependent inelastic deforma- 
tions under particular loading conditions with consequently different response 
characteristics • 

There are more than ten unified constitutive theories in the literature. 
A few of them were proposed in the last five years. These sets of constitu- 
tive equations have some common properties and some essential differences 
which have been reviewed recently by Walker [6] . Since then, there have been 
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more advances in the development of the unified theories. The purpose of this 
survey is to update Walker's previous work by reviewing the state-of-the-art 
and the numerical integration techniques for these unified theories. 

This survey also serves to identify areas for further model developments. 

The unified theories which are reviewed in this survey include those of 
Walker [6], Bodner and Partom [7-9], Miller [10-12], Krieg. Swearengen and 
Rhode [13], Chaboche [14], Robinson [15], Hart and co-workers [16,17], 

Stouffer and Bodner [18,19], Lee and Zaverl [20], Ghosh [21], and Kagawa 
and Asada's modification of Miller's model [22], 
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2.0 GENERAL CHARACTERISTICS OF UNIFIED CONSTITUTIVE EQUATIONS 
FOR ELASTIC-VISCOPLASTIC MATERIALS 


Constitutive equations for elastic-viscoplastic material could be formu- 
lated either with or without the use of a yield criterion. A basic assumption 
for this class of constitutive theories is that in the range where inelastic 
strains are present, the total strain rate e can be divided into elastic and 
inelastic components which are both nonzero, i.e. 


e . . 
ij 



( 1 ) 


This equation is applicable for the small strain case and a similar decompo- 
sition is assumed to hold for the deformation rates in the case of large 
strains. These are equivalent to strain rates if the strains are small. 

For elastic response the stresses are directly related to the deforma- 
tion gradients with no memory effects. The elastic response is fully re- 
coverable both thermodynamically and geometrically. In the case of small 
strains, as considered here, the elastic strain rate is given by the time 
derivative of Hooke's Law. 

In the literature, see e.g. [23,24], there are alternative definitions 
for e.P, the inelastic strain. For our purposes here, and in the context 

p 

of a "unified" theory, the inelastic strain rate will be considered to 

include all strains that are not elastic, i.e. the difference between the 

0 

total strain rate, and the elastic strain-rate e.., (Eq. 1), Thus, the 

ij ij 

expression '’unified” applied to the theories reviewed is taken to mean that 
all aspects of inelastic behavior such as plastic flow, creep, and stress 
relaxation characteristics for different loading histories. This broad de- 
finition of "unified” would admit theories with or without a yield criterion 

and with alternative specifications for e.?. Alternative definitions of the 

1J 

inelastic strain increment that depend on a residual strain after loading and 
unloading from a stress increment are discussed by Lee [23] and Drucker [24]. 
Particular definitions of this type lead to the "normality” condition associated 
with strong material stability. With both loading and unloading increments time 
and temperature dependent, determination of an equilibrium residual plastic 
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strain component is difficult. Alternate uniqueness and stability criteria 
for unified theories with internal variables are discussed in Section 2.5. 

It is noted that the precise definition of inelastic strain rate in a con- 
stitutive theory is given by the equations themselves. 

Constitutive theories which are formulated without the use of a yield 
criterion include that of Bodner and his associates [7-9], Walker [6], Miller 
[10-]2], Krieg, Swearengen and Rhode [13]. Since these models do not con- 
tain a completely elastic regime, the function that describes the inelastic 
strain rate should have the property that the inelastic strain rate be very 
small for low stress levels. 


For theories with a yield criterion, e . . is identically zero until an 

invariant function of the stress reaches a prescribed value; the function, by 

definition, is independent of strain rate. For stresses at or exceeding the 

yield value, Equation (1) applies and e.? and the stresses a., are function- 

# p 1 3 

ally related. The fully elastic state, i.e. = 0, would apply only for 

stress states less than the rate independent yield value, and loading and un- 
loading paths above that are controlled by the loading conditions through the 
constitutive equations. Theories of this type have been developed by Perzyna 
[25] for the case of isotropic hardening and by Chaboche [14], Robinson [15] 
and Lee and Zaverl [20] for the case of both isotropic and directional harden- 
ing. In these theories, loading above the elastic limit value and unloading 
to it would generally not be fully elastic. 


All the unified models are formulated on the basis of internal or struc- 
tural variables which depend on the loading history. The essential features 
of these unified theories are: (1) a flow law which functional form depends 

on the method of treatment of directional (kinematic) hardening, (2) a kinetic 
equation which is the temperature dependent functional relationship between 
the strain rate and stress invariants and includes internal variables, and 
(3) a set of evolution equations for describing the growth of the internal 
variables. Here, the internal variables are used to represent the current 
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resistance to inelastic flow of the deformed solid. Two deforming solids with 
identical values of their internal variables would have identical inelastic 
responses under the same imposed stress state. Both the choice and the number 
of internal variables vary with the unified models. The number of independent 
internal variables that have been suggested ranges from one to as many as six 
reflecting the relative complexity of the model. Most of the unified models, 
however, use two internal variables or one variable with two components: one 

to represent isotropic hardening and another to represent directional (kine- 
matic) hardening. In most models, the isotropic hardening variable is repre- 
sented by a scalar quantity, either the drag stress (K) or the yield stress 
(Y), while directional hardening is represented by a second order tensor 
or a scalar function of such a tensor. 

2.1 Basic Flow Laws 


Four basic forms of the inelastic flow law have been identified. Plas- 
tic incompressibility is always assumed and these flow laws are: 


(1) 

c .? 
ij 


X-S. . 
1 IJ 


# 

e P = 
8 kk 

0 

(2a) 

(2) 

e .? 

ij 

= 

X„£. . 

2 ij 

= 

V s i 

. - o..) 

J iJ 

■ - 0 

(2b) 

(3) 

e , ? 
ij 

= 

X. ... 

ljkl 

s k i 

, e 

P 

iikl 

i* -0 

ljkk 

(2c) 

(4) 

• P 


df 


• P 



(2d) 

e . . 
ij 


do. . 
ij 


* e kk 

= 0 


where S . . , a . . , and 
iJ 

respectively. The 

. are the deviatoric, direct and effective 
ij 

tensor fl. . represents the ff equil ibrium stress" 

stresses, 
which has 

also been 

referred 

to as 

the 

"back 

. tl 

stress 

and the u rest stress”. 

The para- 


meter f is a yield function or a flow potential. It should be noted that the 
first three laws can be considered or can be derived from Equation (2d) if 
they are associated with a flow potential. 

Equation (2a) is the Prandtl-Reuss flow law associated with the von 
Mises yield criterion. However, it can be considered as a basic material 
equation in its own right independently of a yield condition. As such, this 
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equation is usually taken to be applicable for proportional loading conditions 
for which isotropic hardening would be appropriate. The equation states that 
the material response (i.e., the plastic strain rate) to stress is isotropic 
even though could be stress history dependent. Since stress is direc- 
tional, X j could have a directional character within the context of incremen- 
tal isotropic response and thereby account for induced directional hardening 
effects. This situation would be analogous to finite strain elasticity in 
which the coefficients are stress history dependent. 


Equation (2b) is the flow law obtained by introducing the kinematic 
hardening variable of Prager into the classical plasticity formulation to ac- 
count for directional hardening (the Bauschinger effect) [26]. In this con- 
text, the term would represent the new origin of a translating von Mises 
yield surface in deviatoric stress space, and Equation (2b) would be the asso- 
ciated flow rule. As before. Equation (2b) can be taken to be a basic mate- 
rial equation in a formulation without a yield criterion and the ’’equilibrium 
stress tensor A^ is generally intended to serve the following functions: 

(a) to account for directional hardening (the multi-dimensional 
Bauschinger effect), and for the non-coaxiality of and under non- 
proportional loading histories (Figure 1); 


(b) to account for reversed plastic straining effects, e.g. reversed 

creep, relaxation through zero stress, when the effective stress I is neca- 
* . 1 J 

tive; 


(c) for theories without a fully elastic range (i.e., a yield criteri- 
on), to account for low plastic straining within a given range. 

Equation (2b) has been used by many investigators as the basic flow rule 
and appears to be particularly useful in representing directional hardening 
effects. In Equation (2b), X ^ is a scalar function of the isotropic hardening 
variable K and the tensor A... The evolution equation for A. . is generally 
dependent on a hardening term in the direction of and on so called "dynam- 

ic" and "static" (or thermal) recovery terms in the direction of Q . As 

ij 

such, is a deviatoric quantity, i.e. traceless. 

Although attempts have been made to interpret 11.. physically in terms of 
residual stresses, this turns out to be difficult but is not essential from 
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DEVIATORIC stress space 



FIGURE 1. GRAPHICAL REPRESENTATION OF THE BASIC FLOW LAWS USED IN 
THE UNIFIED CONSTITUTIVE THEORIES. For theories based 
on an equilibrium stress, the inelastic strain rate 
vector e!?. is coaxial with the effective stress E-jj and 
normal to J the flow potential f if such a concept is 
used. For theories which do not include an equilibrium 
stress, is coaxial with the deviatoric stress Sjj 
for both isotropic and incremental ly isotropic cases but 
is noncoaxial with S-jj for generalized anisotropic cases. 
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the viewpoint of its possible utility in constitutive equations. The term 
"equilibrium stress" is usually applied to since it would correspond to 
the asymptotic stress state under relaxation conditions, but such a stress 
state is difficult to determine by direct experimentation. A possible criti- 
cism of Equation (2b) is that the "equilibrium stress" is a load history de- 
pendent material property which is subtracted from the applied stress which is 
a kinetic quantity. It is noted that difficulties are experienced in general- 
izing the flow law (2b) to the case of large strains because different trans- 
formation rules are required for and for , e.g. recent papers by E. H. 

Lee [27]. 


For cycling under proportional stress states. Equation (2a) with a 
stress history dependent coefficient can be shown to be equivalent to Equation 
(2b). The real differences in these equations would show up for non-propor- 
tional loading histories. 

Equation (2c) is the general anisotropic form of the Prandtl-Reuss flow 
law which can be rewritten in a 6D stress and strain rate space to take the 
form. 


E P 

a 


A . T n 

ap 0 



(3) 


where E and T n are related to the usual plastic strain rates and stresses in 
o p 

a simple manner, see [18,19], and A D is the 6x6 matrix of coefficients. If 

ap 

the material is initially isotropic and the law for plastically induced direc- 
tional hardening does not lead to off diagonal terms, then is initially 

and remains diagonal. Under these conditions. Equation (2c) is equivalent to 
Equation (2b) since 6 material constants determine the anisotropic flow be- 
havior (fiL ^ i n Equation (2b) has 6 components) • All the flow equations. Equa- 
tions (2a,b,c), would be equivalent for the case of proportional loading, in- 
cluding cyclic conditions. 


Although Equation (2c) can be used for materials that are initially non- 
isotropic with regard to inelastic response, it is obviously a complicated 
equation to manage. There has been relatively little experience in using 
Equation (2c) for such cases. 
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For constitutive theories with a flow potential, both the flow law and 
the growth law of the kinematic hardening variable 0^ are derivable from a 
single flow potential. Equation (2d). A basic form of such a flow potential 
is [13,28]. 


f = f (o^. Q^) « F(J 2 ) + G(T 2 ) 


(4) 


where J 2 = 1/2 (o^ - 0^) (o^ - ) 


T 2 - 1/2 °ij °ij 


If the function F(J 2 ) taken as [15] 


n+1 


F(J 2 ) - 


K 


P (n+1) 


(5) 


j2 

where F = — ^ - 1 is the von Mises yield function, n and p are material para- 
meters, and recalling that K is the isotropic hardening variable. The associ- 
ated flow law, from Equation (2d) then becomes 

n-1 

= F ^ ” Qjj ) for inelastic loading (6a) 

and a.? = 0 for elastic unloading (6b) 


The conditions for inelastic loading and elastic unloading have been identi- 
fled in [15]* It can be easily seen that (6a) and (2b) are equivalent* In 
both cases* the direction of the inelastic strain rate vector is coaxial with 
the current effective stress vector (see Figure 1). Robinson [15] has recog- 
nized the equivalence and pointed out that the existence of a flow potential 
is not strictly necessary for the development of the flow law represented by 
(6a); the flow law can simply be stated without reference to their derivabil- 
f r °m a potential function* as done in models which do not include a yield 
surface or flow potential* 


2 .2 Kinetic Equations 

The flow laws* Equations (2a) and (2b) can be squared to give respec— 
tively , 

l l - (7*) 


A- 1 0 


(7b) 


s - 'V'v 1 ' 2 


where D. is the second invariant of the plastic strain rate, D P = (1/2) s P 
. p ' 2 ij 
s^j* and and are the second invariants of the deviatoric stress and 

effective deviatoric stress* respectively* 


- (1/2) Sy Sy 

J 2 * (1/2 » <S ij - V‘ S ij - n ij ) 


(8a) 

(8b) 


Fundamental to all ’’unified" viscoplastic formulations based on flow 

laws of the forms listed in Equations (2) is that inelastic deformations are 

governed by a functional relation between D y and J (or J ) that could in- 

X z z 

volve load history dependent variables. These variables are intended to re- 
present properties of the inelastic state with respect to resistance to plas- 
tic flow, e.g. hardening and damage. Some functions that have been suggested 
are the following. 


(a) 

D 2 P ’ V" 



(9a) 

(b) 

d 2 P - D o exp 



(9b) 

(c) 

D/ - D o [,UJ,(I>"]“ 


(9c) 

where 

X = 3J 2 /K 2 

* or 3 J 

* 2 
2 /K 


and n, at. 

and D are constants. The inelastic strain rate components 

can then 


be obtained as a function of the stress by the use of Equation (2a) or (2b) 

and one of Equations (9). Expression (9b) has some advantage over (9a or c) 

in theories without a yield criterion in that the value of is almost zero 

for some extended range of J. regardless of the value of n. In (9b)* D is 

i o 

the limiting value of the inelastic strain rate in shear; (9a) and (9b) do not 

contain such a limit. These differences between the kinetic equations are 

illustrated in a normalized plot of log (D-^/D ) vs X in Figure 2 for the case 

z o 

of n = 3 and m = 1.0. Experimental support for Equation (9b) is shown in 
Figure 3 which indicates that a limiting inelastic strain rate of 2x10^ sec ^ 
occurs in 1100-aluminum [29]. On the other hand* the power law expression has 
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KINETIC RELATIONS 



FIGURE 2. FUNCTIONAL BEHAVIOR OF THE KINETIC EQUATIONS USED IN THE 
UNIFIED CONSTITUTIVE THEORIES. The exponential formula- 
tion in Bodner-Partom's theory is seen to give a limit- 
ing inelastic strain rate of Dp. 
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FIGURE 3. STRAIN RATE AS A FUNCTION OF STRESS FOR 1100 ALUMINUM 
INDICATING A LIMITING STRAIN RATE OF 2 x 105 SEC-1, 
FROM [27] 




been found to overpredict values of stress in a constant strain rate tensile 

—2 —1 

test when the strain rates are greater than 10 sec [13,30]. 

In all the preceding equations (9a,b,c) the exponent n influences the 
slope of the D^, relation and therefore has the major influence on strain 
rate sensitivity. That parameter also affects the overall level of stress- 
strain curves although the level also depends on the hardening parameter E. 

Temperature (T) dependence of plastic flow is a first order phenomenon 
comparable to strain rate sensitivity and should appear directly in the ki- 
netic equation. In the case of Equations (9a, b), this can be achieved by 
taking the exponent n to be a function of T, e.g. n = c/kT (k is Boltzmann's 
constant and c a material constant) which leads to strong temperature depend- 
ence of the stress parameter X=3J 2 /K 2 (or 3J 2 /K 2 ) . Numerical results for this 
dependence are shown in Figure 4 for both the power law and exponential kinet- 
ic equations at different non-dimensional ized strain rates. These trends ap- 
pear to be consistent with experimental results shown in Figures 5 and 6 which 
are plotted in a similar manner. 

The method of including temperature dependence in Equations (9) can be 
derived from an activation energy formulation. Table I lists temperature- 
dependent kinetic equations based on four different functional expressions for 
the activation energy. Some of the consequences of the various relations are 
discussed in [29] . 

Another procedure for including temperature dependence in the kinetic 
equations is to multiply the stress function, the right hand side of equa- 
tions (9) by a temperature function. The temperature factor can again be 
motivated by thermal activation considerations and the Arrhenius expression 
seems to be the reasonable function to use at the higher temperatures. This 
is the approach taken in [10]. Additional discussion of temperature effects 
in the constitutive equations is given in a subsequent section. 

2.3 Evolutionary Equations for Internal Variables 

An important feature of the unified approach is the use of a set of dis- 
crete evolutionary equations to describe the change in hardening behavior of 
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FIGURE 4. FUNCTIONAL BEHAVIOR OF TEMPERATURE-DEPENDENT KINETIC 
EQUATIONS UTILIZED IN BODNER-PARTOM AND WALKER 
THEORIES 
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TABLE I 


FIVE FORMS OF TEMPERATURE-DEPENDENT KINETIC EQUATIONS WITH 
THE CORRESPONDING ACTIVATION ENERGY FUNCTION 


Activation Energy Temperature-Dependent Kinetic Equations 


Ml = C£n 


3J, 


D 2 * D 0 


r3J 2 i c/kT 


L K 


= H o - Vg(J 2 ) 


D 2 = D 0 ex P 


H 0 " Vg( 


kT 


AH = 


H*r 

3Jo 


D 2 = D 0 ex P 


H* i 

K 

kT ' 

3J„ 


AH = kT 


[~K 2 1 

C/kT 


/ ? \C/kT1 

3J 2 


D 2 = Dgexp 

1 1 

^coil 
° rf 

ro 1 ^ 

i 


AH = Q 



where C, D 0 , H*, H 0 , Q, m, and n are constants; 
volume; and k is the Boltzmann's constant. 


V is the activation 
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material s undergoing inelastic deformation. The same set of internal vari- 
ables and evolutionary equations is used to govern all aspects of inelastic 
deformation including plastic flow, creep, and stress relaxation. The general 
framework of the evolutionary equations of internal variables is based on the 
now well-accepted Bailey-Orowan theory [31-32] which theorizes inelastic de- 
formation to occur under the actions of two simultaneously competing mech- 
anisms, a hardening process proceeding with accumulated deformation and a re- 
covery or softening process proceeding with time. The evolution rate of an 
internal variable is then the difference between the hardening rate and the 
recovery rate as given by 

ij - w *1 - w T) <10> 

where Xj is the evolution rate of the internal variable X^ and hj and r x are 
the hardening and the thermal recovery functions, respectively. and x ± are 

functions of X., temperature (T) . and the hardening measure (Mj) is either 
J p or I p depending on the model. A summary of the hardening measures uti- 
lized in the various unified theories is shown in Table II. 

(1 .) Isotropic Hardening 

The quantity K in Equation (9) is usually interpreted as the iso- 
tropic hardening internal variable and is often referred to as the drag 
stress. Evolutionary equations for the isotropic hardening parameter general- 
ly follow the hardening/ recovery format shown in Equation (10). A comparison 
of these hardening and recovery functions in various unified theories is shown 
in Table III. The rate of isotropic hardening is usually given by a function 
of the hardening variable X, which may saturate to a limiting value, shown as 
K in Table III, multiplied by a measure of the hardening rate. Both the in- 
elastic work rate and the effective inelastic strain rate have been proposed 
as the scaler hardening measure. At present, it is not obvious which harden- 
ing measure leads to better overall agreement with experiments. A recent pro- 
posal is to take the measure to be a functional of the strain rate history 
[33], On the other hand, the rate of softening or recovery is often taken to 
be a power function of K and a temperature-dependent constant K q which value 
represents the reference state for that particular temperature. 
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THE CHOICE OF INTERNAL VARIABLES AND MEASURES OF 
HARDENING IN SELECTED UNIFIED CONSTITUTIVE THEORIES 
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MILLER 



TABLE III 


THE SPECIFIC FORMS OF ISOTROPIC HARDENING AND STATIC THERMAL RECOVERY 
FUNCTIONS USED IN THE SELECTED UNIFIED CONSTITUTIVE THEORIES 

K = h-j (K)M-j - r ] (T,K) 

where M-j = e; e = V J e ij ^ij 

• • 

or M-| = Wp (Bodner-Partom's Theory) 

Static Thermal 

Model Hardening Function, h](K) Recovery Function, n(T,K) 


Bodner-Partom 

- K) 

c 2 (K-K 0 ) p 

Wal ker 

- K) 

- 

Krieg et al 

C 1 

c 2 (K-K 0 ) p 

Robinson 

C 1 

- 

Chaboche 

C-| (K-j - K) + f 1 (S, e, a..) 

- 

Lee and Zaverl 

C-| (K* - K)/ST Z 

- 

Hart 

- 

- 

Ghosh 

C, 

c 2 (k-k 0 ) p 

Miller 

C, [K, - C 4 (sinh _1 C 3 |e| ) q K m ] 

C 2 [sinh C 3 K m ] 


* 

where C-| , C 2 , C^, C^, C^, m, p, q, K Q , and K-j are material constants; K-| 
is the saturated value of K; K* is governed by an evolutionary equation 
which is function of c and J 2> 
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This recovery model, sometimes credited to Friedel [34], theorizes that 

recovery occurs only when the current internal state exceeds the reference 
state • 

Directional or Kinematic Hardening 

Probably the main difference in the various unified theories is 
the treatment of directional or kinematic hardening. Differences exist not 
only in the choice of the flow law but also in the evolutionary equations. 

The general framework of these evolutionary equations follows the hardening/ 
recovery formulation represented in Equation (10) with indexes to indicate 
the directions of hardening and recovery. 

°ij = h 2 ( °ij )A ij " d(Q ij’ T) *ij ~ r 2 (Q ij'' r)V ij + e <° ij * T) * W ij Ul) 

where h^, d, and r 2 are the hardening, "dynamic recovery", and static thermal 
recovery functions, respectively. 0 represents hardening and/or recovery as- 
sociated with the rate of temperature change. M... N, , V,. and W. are the 

ij ij ij ij 

directional indexes of h 2 > d, r 2 , and 0, respectively. The main differences 
among the various theories, as summarized in Table IV, are in the choices of 
the directional index and the hardening and recovery functions. 

As indicated in Table IV, unified models based on the equilibrium 
stress utilize the inelastic strain rate as the directional index for harden- 
ing and contain a "dynamic recovery" term in the hardening function. The pro- 
posed hardening rule is thus similar to the Prager rule [26] in conventional 
plasticity which requires the translation of a yield surface to occur in the 
direction of the plastic strain increment. On the other hand, the evolution- 
ary equation proposed in conjunction with Equation (2a) is based on the direct 
stress as the index for directional hardening [8,9]. This formulation avoids 
the cross-softening effect associated with inelastic strain rate as the index 
and the theory is more compatible with Ziegler's modification [35] of the 
Prager hardening rule. 

The directional index for "dynamic recovery" is generally in the 
opposite direction of the directional hardening variable . The "dynamic 
recovery term is treated in [7-9] as a saturation term in the direction of 
the direct stress but the index has recently been modified to be in the direc- 
tion of - flij also [36] . 
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THE SPECIFIC FUNCTIONS OF ANISOTROPIC HARDENING, DYNAMIC RECOVERY, STATIC THERMAL 
RECOVERY, AND THE TEMPERATURE RATE TERM IN SELECTED UNIFIED CONSTITUTIVE THEORIES 
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are the saturated values of ft.,; fl!, are governed by evolutionary equations which are functions 



The unit vector which represents the direction cosines of the 
directional hardening variable is usually taken to be the directional index 
for static thermal recovery. Recovery always occurs in the opposite direction 
of the unit vector and tends to reduce the magnitude of the directional (kine- 
matic) hardening variable. Most unified theories utilize Friedel's recovery 
model and take zero magnitude of as the reference state. Table IV shows 

that a temperature rate term is also included in the theories of Walker and 
Chaboche • In principle, similar terms could be added to the other theories. 

The various hardening terms in the evolution equation of the 
directional hardening variable have profound effects on the shape of the cy- 
clic hysteresis loops. Unified models which contain at least two hardening 
terms to describe the evolution of the equilibrium stress with strain history 
include that of Walker 16 ], Hart and co-workers [17], and Miller [12]. The 
contribution to the equilibrium stress from each of the four hardening terms 
used in the Walker theory may be seen in Figure 7. The term describes the 

difference between the tensile and compressive behavior of the equilibrium 

stress variable, and this difference gives a stress-strain response which 
differs in tension and compression. The linear term n^e.^ allows the asymp- 
totic stress-strain response to exhibit a linear behavior at large strain val- 
ues. The integral term containing allows the equilibrium stress to grow and 
saturate rapidly when cycling under small strain amplitudes, while the inte- 
gral term containing Ug gives a more gradual growth of the internal variable. 
These two integral terms allow the cyclic hysteresis behavior to be reproduced 

for both small and large strain amplitude loading conditions. 

2.4 Temperature Dependence of Inelastic Flow 

Temperature effects are taken into account in most unified models by 
assuming the material constants in the kinetic equation and the evolutionary 
equations of the internal variables to be temperature dependent. Lindholm and 
Bodner [29] recently examined this issue and made an attempt to model tempera- 
ture dependence of inelastic flow by incorporating a reaction rate theory into 
Orowan's kinetic relation. By assuming appropriate functions for the activa- 
tion energy, Lindholm and Bodner were able to obtain four different formula- 
tions of temperature-dependent kinetic relations. These kinetics relations 
and the corresponding activation energy are rewritten in terms of deviatoric 
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invariants, summarized in Table I. Figure 4 demonstrates the similarities 
between the functional behavior of the first and fourth formulations which 
correspond to the kinetic relation utilized in the Walker and Bodner models, 
respectively. Experimental evidence in Figures 5 and 6 for aluminum [29] and 
iron [37], respectively, tend to indicate both formulations are adequate and 
justified. 

Another approach for incorporating temperature dependence into a unified 
theory is taken by Miller [10], Temperature dependence of inelastic flow is 
modeled on the basis of creep^ phenomenology and an Arrhenius expression is 
modified to reflect the temperature dependence of the activation energy for 
diffusion (or creep). To cover a wide range of temperature, both lattice dif- 
fusion and pipe diffusion have been considered. The resulting Arrhenius ex- 
pression which is also shown in Table I is included in both the kinetic rela- 
tion and the recovery term of the evolutionary equations. With the exception 
of the elastic constants, material constants in Miller's model are independent 
of temperature. There are, however, indications that the model would yield 
better results for thermomechanical loading if some of the material constants 
are allowed to change with temperature [6] . 

Walker and his co-workers at Pratt S Whitney Aircraft have compared the 
predictive capabilities of three unified models for applications under thermo- 
mechanical cycling conditions. Their experience with the unified models to 
date indicates that all the material constants in the formulations would de- 
pend on temperature and must be evaluated at a number of base temperatures. 
Suitable interpolation methods must be developed for evaluating the values of 
these material constants at other temperatures. The predicted results would 
therefore depend on the interpolation method [6]. 

The temperature dependence of the internal variables is also important. 
Comparisons of experimental and theoretical results on thermomechancial con- 
stitutive behavior of Hastelloy-X in Reference 6 indicate that the shape of 
the predicted thermomechanical loop depends critically on the growth of the 
equilibrium stress with temperature. Most thermomechanical loops are neces- 
sarily of small strain amplitude. Hence, at the points of strain reversal, 
the equilibrium stress is in the initial growth phase and is changing rapidly 
with strain. Depending on the rapidity with which the equilibrium stress 
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grows with inelastic strain after a strain reversal, the equilibrium stress 
may be in the process of almost saturating at the next strain reversal. This 
{in change the predicted shape of the thermomechanical stress - strain hyster- 
esis loop quite dramatically. 

2.5 Uniqueness and Stability Criteria 

All the proposed kinetic and evolutionary equations of the internal 
variables are based on a combination of phenomenology and physical mechanisms. 
There are certain continuum properties which are required of these constitu- 
tive relationships in order to ensure that the resulting boundary value prob- 
lems are well-posed. The required properties are the uniqueness and the 
stability of the solution. For stability, unified theories with internal 
variables must, according to Ponter [38], obey the following inequality: 


do. . de . 


ij ij 

where do.., de ® 


- dX i dX £ >, 0 


( 12 ) 


U 6 -j» dX i» and dX i represent incremental changes in stress, in- 
elastic strain rate, the current value and the evolution rate of the internal 
variables. The inequality admits classical plastic flow, creep, and stress 
relaxation behavior. It also admits recovery phenomena involving negative 
inelastic work provided that the corresponding changes in the internal vari- 
ables are sufficiently large to make the inequality in Equation (12) remain 
valid. The basic requirement of Equation (12) is that the dissipation rate 
must be nonnegative. 

For a constant internal state, a small change in results in a cor- 
responding change in so that [38] 


d<T ij d *ij >/0 ' *i = ° U3> 

The inelastic work inequality which is identical to Drucker's postulate [39] 
in classical plasticity requires that for a stable material flow the plastic 
work done must be nonnegative. For proportional loading the kinetic equations 
represented in Equation (9a) to Equation (9c) all yield convex flow poten- 
tials to which the inelastic strain rate vectors are normal . The consequence 
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is that the inelastic work is always positive, and unified theories based on 
Equation (9a) to (9c) obey the inelastic work inequality. 

For loading involving constant as in creep, the inequality in Equa- 
tion (12) requires that a small change in dX^ produces a change in the evolu- 
tion rate such that 


- dX £ dX £ > y 0 , do^ = 0 (14) 

It is clear that both Equation (13) and (14) would be obeyed by any constitu- 
tive model if these two conditions are satisfied: 


( 1 ) 

( 2 ) 


An increase in strain rate (de^ >/ 0) results in an 
increase in flow stress (da^j ^ 0) or vice versa; and 

An increase in the value of an internal variable 
(dX^ *>/ 0) results in a decrease in the evolution 
rate (dX^ ^ 0) or vice versa. 


Some of the ramifications of Drucker's inequality in classical rate 
independent plasticity are that the stress-strain curve must have a positive 
slope for stable flow and strain must be a single-valued function of stress 
for uniqueness. Ponter's inequality exerts similar ramifications on visco- 
plastic unified theories. For uniqueness, it appears that the inelastic 
strain rate must be a single— valued function of stress and internal variables. 

To satisfy the requirement for stable flow, Equation (12) dictates that stress- 
strain curves at constant strain rate must have positive slope but must decrease 
with increasing strain. On the other hand, stress-strain curves at constant 
value of plastic strain or plastic work must have positive slope, but the slope 
may either increase or decrease with increasing strain rate [38]. 


Most (if not all) of the unified theories listed in Table II satisfy the 
Ponter inequalities and meet the uniqueness and stability requirements. The 
stability requirement is, however, not essential for constitutive theory de- 
velopments. Unified theories admit unstable inelastic flow and are generally 
modeled by including softening mechanisms such as thermal softening and con- 
tinuum damage in the evolution and/or kinetic equations. Negative strain rate 
sensitivity and strain (or work) softening phenomena have also been modeled 
[ 11 , 12 ]. 
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2.6 Constitutive Behavior Under Nonproportional Loading 

Cyclic constitutive behavior of materials subjected to nonproportional 
loading conditions is still poorly understood. Recent studies on Hastelloy-X 
[36], copper [40], 1% Cr-Mo-V steel [41], and 316 stainless steel [42] all 
indicate that these materials exhibit considerably more cyclic hardening when 
tested under nonproportional paths of combined tension and torsion than under 
proportional paths of tension or torsion only. A comparison of the hardening 
behavior of Hastelloy-X alloy subjected to tension, torsion, and combined, 
nonproportional tension-torsion loading is shown in Figure 8 [36], In this 
plot of effective stress vs cumulative effective plastic strain, where the 
effective stress and strain are computed based on the von Mises criterion, and 
the cumulative effective plastic strain is number of cycles times the effec- 
tive plastic strain range per cycle, an apparent increase of hardening is ob- 
served in the case for 90° out-of-phase loading. This apparent hardening can 
arise from three different sources: (1) nonproportional loading induces addi- 

tional hardening; (2) the von Mises criterion is not appropriate for this ma- 
terial; and (3) testing was not done at the same "effective” strain levels. 
Efforts to delineate the relative contribution of these factors indicate the 
first factor to be important. As a result, all of the constitutive models 
need to be modified to take into account the hardening behavior due to out-of- 
phase loading. 

Another question concerning the nonproportional, multiaxial deformation 
is whether or not the inelastic strain rate vector is coaxial with the devi- 
atoric stress vector under nonproportional loading paths. The works of Ohashi 
and co-workers [43-45], and Meguid [46] indicate that the deviatoric stress 
and the inelastic strain rate vectors become noncoaxial just after an abrupt 
change of deformation path but become coaxial again as inelastic deformation 
develops along the subsequent path. Current work by the authors on 90° out- 
of-phase cyclic loading of Hastelloy-X under combined tension and torsion in- 
dicates that the inelastic strain rate vector leads the deviatoric stress 
vector by an angle which ranges from 10 to 40 degrees [36] . These findings 
tend to suggest that, as a result of the apparent noncoaxiality, the isotropic 
Prandtl-Reuss flow laws might not be adequate for nonproportional loading 
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Effective Stress, MPa 



Cumulative Effective Plastic Strain, % 


FIGURE 8. A COMPARISON OF THE EFFECTIVE STRESS-CUMULATIVE EFFECTIVE 
STRAIN BEHAVIOR OF HASTELLOY-X UNDER TENSION, TORSION, 

AND COMBINED, NONPROPORTIONAL TENSION AND TORSION 
LOADINGS 
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paths. On the other hand, the noncoaxiality between inelastic strain rate and 
dewiatoric stress can be accounted for by using either a generalized anisotro- 
pic Prandtl-Reuss flow law or a flow law based on an "equilibrium stress'.' 
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3.0 REVIEW OF NUMERICAL METHODS FOR INTEGRATING 
UNIFIED CONSTITUTIVE EQUATIONS 


3*1 Numerical Integration Techniques 

A literature review indicated that no systematic studies which compare 
the integrabil ity of unified constitutive equations have been performed. The 
evaluations that have been done support the existing consensus that this fam- 
ily of differential equations can be characterized as mathematically "stiff," 
That is, in these equations, dependent variables are susceptible to large 
changes from small increments of the independent variables or from small time 
steps. This "stiff" behavior occurs usually with the onset of a significant 
amount of inelastic strain in the loading cycle and is due to the generally 
nonlinear nature of the functional forms that are employed in the kinetic 
equations of these theories. 

Probably of greater importance than differences in stiffness of indi- 
vidual constitutive theories are the methods used to integrate them. A sys- 
tematic comparison of a variety of approaches has been reported by Kumar, 
Morjaria, and Mukherjee [47]. Specifically, they compared simple Euler for- 
ward method, two-step Adam's method, predictor-corrector method, and the Gear 
method. This study concluded that for the constitutive theory of Hart, a re- 
latively simple Euler integration method, together with a time step control 
strategy, was optimal when compared with the more sophisticated methods. Sim- 
ilar conclusions were obtained by Krieg [48] who rearranged various unified 
constitutive theories into a skeletal model to illustrate the numerical dif- 
ficulties in the integration of constitutive equations and to discuss the 
viability of various integrators. 

Tanaka and Miller recently developed a noniterative, self-correcting 
solution (NONSS) method for integrating stiff time— dependent constitutive 
equations [49]. The NONSS method is basically an extension of the a— method 
which is used in creep and heat transfer analysis [50]. In this approach, the 
stress, inelastic rate, and evolutionary rate of the internal variables at 
t+oAt are expressed as: 


cr(t+aAt ) 


(l-a)a(t) + aa(t+At) 


i(t+oAt) 


(l-a)e(t) + ae(t+At) 


(15) 


i^t+aAt) * (l-o)i 1 (t) + aX^t+At) 

with 0 ^ a ^ 1. The method reduces to the Euler forward difference method 
when a = 0 and is an implicit technique when a > 0, because the unknown quan- 
tities cr( t+At) , e(t+At), and X^t+At) enter into the solution. Implicit quan- 
tities are removed in the NONSS method by Taylor expansions of c, and 
X^. The method is unconditionally stable for a!^ 1/2 and noniterative, but 
requires setting up of a Jacobian matrix and solving a set of linear equations 
at each time-step. Accuracy is maintained through self-adaptive time control 
and by correcting previous errors at the current step. This method has been 
used in one-element applications only. The applicability of this approach to 
finite-element analysis remains to be seen. 

At Pratt and Whitney Aircraft and at the United Technologies Research 
Center, work has been done using the Euler single step approach usually with- 
out automatic time step control, but rather by determining an optimum step 
size for each problem experimentally. Efficiency obtained by using this ap- 
proach has been acceptable and has shown considerable improvement over more 
sophisticated approaches such as higher order Runge-Kutta methods. 

One way of incorporating a self-adaptive time control technique into the 
integration procedure is by expressing the inelastic strain at time t+At in a 
Taylor series and ensuring that the higher order terms are small in comparison 
to the first order term. For example, the inelastic strain at time t+At may 
be written in the Taylor series 

e(t+At) = 6 (t) + e At + 1/2 e (At) 2 +... (16) 

The time step At can then be chosen so that the third term in the series is 
some small fraction* y say, of the second term in the series. In this way one 
obtains 

A+ _ Zx 111 (17) 

At - <17) 

where lei and lei are the magnitude of the maximum components of e and e, 
respectively. 
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Alternatively, implicit methods can be used. An implicit method for use 
with an integral constitutive formulation was outlined in Reference 6 and the 
technique gave good results. 

A summary of these various numerical techniques and their applications 
to several unified theories as well as to Norton's law for integrating a uni- 
axial stress-strain curve to a total strain of 1-2% is shown in Table V. As 

illustrated in Table V, the explicit Euler method is stable when the size of 

-4 

the strain increment is kept below 10 • The size of the strain increment can 

be increased by using an implicit method such as the NONSS or a-method with 
a - 1 (implicit Euler). By restricting the comparison to the explicit methods 
only, it appears that there is no substantial difference between the integra- 
bility of Walker [6,51] and Miller [52] theories nor between these unified 
theories and the classical Norton law [52] . The size of the strain increment 
is, however, somewhat sensitive to the values of model constants which de- 
scribe material strain rate sensitivity. 

3 .2 Integration of Unified Theories in Finite Element Analysis 

The question of which integration method to use for integrating unified 
theories in finite element analysis appears to be code and problem dependent. 
For example, the MARC code solves the finite element equilibrium equations by 
a Picard, or successive substitution, method* The finite element equilibrium 
equations may be written in the form 

£/ B T <i(u)dV = P (18) 

P being the applied load sector and o(u) being the corresponding stresses due 
to the nodal displacement u. 

The MARC program solves the incremental form of equation, viz. 

^z/ v B T DBdV^ Au i+1 = AP - Z / y B T A?(Au i )dV. (19) 

where is the inelastic stress increment in the relation 

Ao = DAe + A^ , (20) 
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COMPARISONS OF THE INTEGRATABILITY OF VARIOUS CONSTITUTIVE MODELS 
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Average strain increment per step = total strain/number of time steps. 
























are 


and Au^ + ^ are the nodal displacement increments in the (i+1) the interation. 
The load AP in Equation (19) represents the current load increment plus the 
out-of-equilibrium load from the previous increment. 


For crack problems the preceding interation method can converge very 
slowly and very small finite element load steps must then he taken. The 
ABAQUS code solves the finite element equilibrium equations directly by a 
Newton-Raphson method, so that Equation (18) is solved in the form 

Bd ' r 


c 1+ i 


= P -if B 


oU^dV 


( 21 ) 


where 


u i+l = U i + C i+1 

and C^ +1 is the ( i+1 ) th interative correction to the total nodal displace- 
ment at the end of the increment. For crack problems small steps must still 
be taken, but now the convergence of the finite element equilibrium equations 
is quadratic, as opposed to the slow linear convergence of the successive 
substitution method. 

To obtain the quadratic convergence it is necessary to evaluate the 
Jacobian matrix da/de. Since subincremental procedures are usually used to 
integrate the unified models over the finite element increment, the Jacobian 
matrix cannot be obtained analytically. It may be determined by perturbing 
each component of the strain increment tensor and obtaining Bat Be numerically. 
It is then necessary to use the same number of subincrements in each of the 
perturbations so that do/da is consistent and does not include contributions 
from changes in the number of subincrements . In this case it is better not to 
use a self adaptive integration procedure. A maximum strain increment of, say 
10 -4 , can be prescribed. The finite element strain increment can then be di- 
vided by 10 4 to obtain the number of subincrements in the step at the parti- 
cular Gauss integration point. If instability is detected during the finite 
element load step, the integration can be repeated using triple the number of 
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subincrements until the stability criterion is met. This number of subincre- 
ments is then used in each integration of the perturbed strain increment 
tensor to obtain the Jacobian matrix do/de. 

This procedure has been implemented in the ABAQUS code to solve crack 
growth problems with a unified model [53] • Typically at the integration point 
closest to the crack tip, as many as fifty subincrements may be required. 
However, since the strain increments at integration points far removed from 
the crack tip are small in comparison with those at the crack tip, only two 
subincrements may be needed at most of the other integration points in the 
structural model. Since small finite element increments are required to solve 
crack problems, and only two or three subincrements are required at most of 
the Gaussian integration points away from the crack tip vicinity, it may ac- 
tually be computationally inefficient to use an implicit integration method 
for the unified models for this class of structural problems. 

3 .3 Comparison of the Integrabil ities of Unified and Classical Theories 

In order to better understand the problem of efficient computing algo- 
rithms for utilizing unified constitutive models, a summary of computing times 
for typical non-linear finite element analyses has been compiled along with a 
comparison with classical approaches for the same analysis cases. A summary 
of this study is shown in Table VI. This table lists two small structural 
models and some single element models which were performed at Pratt and 
Whitney using the MARC non-linear finite element code. Each of the three 
structural models shown in Table VI was run using two or more approaches for 
modeling the material behavior as indicated. It should be noted that each of 
the variations was run on a single computing system, but that analyses 2 and 
3 were run on a computing system that was approximately 30% slower than for 
the structural model used for analysis number 1, the Burner Liner Crack Anal- 
ysis. Thus, any attempt to make comparisons between structural models may be 
misleading. Aside from this limitation, it can be seen that experience to 
date using the unified model of Walker has involved more computing time than 
the classical approach using a time independent plasticity model and alternate 
load increments using a time dependent creep model. Also shown in Table VI 
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TABLE VI 


RUN TIME COMPARISONS FOR UNIFIED CONSTITUTIVE MODEL APPLICATIONS 


1. Burner Liner Crack Analysis 

Problem Size: 50 elements, 181 nodes, 350 degrees of freedom, 40 load 
increments 


Computing Time 

Method (Seconds) 


a. Classical plasticity/creep 365 

b. Integral form of Walker's theory 322 

c. Walker integral form with 2 subincrements 382 

d. Walker integral form with 5 subincrements 493 


2. Hot Spot Blister Analysis 

Problem Size: 23 elements, 103 nodes, 198 degrees of freedom, 51 or C4 
load increments 


Computing Time 

Method (Seconds) 

a. Classical plasticity/creep 444 

(64 load increments) 

b. Differential form of Walker's model 614 

(51 load increments) 


3. Uniaxial Analysis 


Problem Size: 1 elements, 8 nodes, 12 degrees of freedom, 35 or 48 load 
increments 


Computing Time 

Method ( Second s ) 

a. Elastic Analysis 

(35 load increments) 

b. Classical plasticity/creep 12 

(35 load increments) 

c. Walker's theory, integral form 32 

(48 load increments) 

d. Walker's theory, integral form 48 

(48 load increments with 5 subincrements 

per load step) 
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•re comparisons of computing time for variation of integration method for 
Walker’s constitutive theory. These should not be regarded as a measure of 
efficiency in any absolute sense since no attempt was made to optimize the 
approaches used, but rather are intended to provide some indication of current 
experience for computing times and to provide an informal basis against which 
to judge the computing efficiency of the unified theories. 
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4.0 REVIEW OF SPECIFIC UNIFIED THEORIES 


There are at least 12 viscoplastic constitutive theories which are based 
on the unified approach in the literature. These theories are reviewed indi- 
vidually as follows. 

Other unified theories include that of Lee and Zaverl [20], Cernocky and 
Krempl [30], and Anand [59]. The model of Lee and Zaverl is a generalized 
anisotropic theory based on a yield surface concept. This model and that pro- 
posed by Cernocky and Krempl have recently been reviewed by Walker [6] . 

Anand' s model is in uniaxial form with isotropic hardening only, and its for- 
mulation is very similar to Bodner and Partom's theory. 

4 .1 Robinson's Theory 

Robinson's model [15] is based on a flow potential from which the flow 
law and the growth laws of the internal variables are derived. For stress 
states inside the flow potential, the material is elastic, while for stress 
states on the flow potential, the material response is viscoplastic. Differ- 
ent forms of flow and growth laws are derived for loading and unloading, to- 
gether with inequalities defining boundaries across which the flow and evolu- 
tionary equations change form discontinuously • The model employs an equili- 
brium stress internal variable 0. # and a drag stress internal variable K to 

1J 

represent kinematic and isotropic hardening, respectively. The use of differ- 
ent growth laws allows the equilibrium stress to grow rapidly with inelastic 
strain just after a stress or strain reversal, but to grow more slowly with 
inelastic strain at a large distance (in strain space) from the strain or 
stress reversal. As a result, the model can reproduce rounded, cyclic hyster- 
esis loops for both large and small strain amplitudes without employing two 
hardening terms in the equilibrium stress evolutionary equations. 

The model contains eight temperature-independent and four temperature- 
dependent material constants. The model has been used for predicting material 
behavior under monotonic, cyclic, creep, stress relaxation, and nonisothermal 
cycling loading conditions. [53]. 
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•2 Walker's Functional Theor 


This theory was developed by modifying the constitutive relation for a 
three-parameter viscoelastic solid. Two internal variables were introduced 
into the viscoelastic theory to account for the effects of viscoplasticity. 
The equilibrium stress represents nonlinear kinematic hardening and ac- 
counts for the Bauschinger effect, while the drag stress K represents isotro- 
pic hardening and accounts for cyclic hardening of the material. Both the 
integral and the differential forms of the theory are summarized in Reference 
6 • 


The growth law for the equilibrium stress contains both dynamic recovery 
and static thermal recovery terms. At high strain rates, the thermal recovery 
term becomes insignificant in comparison with the dynamic recovery term, and 
the equilibrium stress becomes independent of strain rate. The equilibrium 


stress expression contains a temperature rate term, which allows the equilib- 

o 

rium stress to change during nonisothermal elastic excursions, and a term 
to account for asymmetric cyclic hysteresis behavior. Static thermal recovery 
terms have been omitted in the growth law for the drag stress. This form has 
been found adequate in the modeling of Hastelloy-X behavior [6], but future 


applications may require the inclusion of static thermal recovery in the drag 


stress evolution law. 


The kinetic relation is formulated in terms of a power law in [6] and in 
terms of an exponential law in [53]. The power law expression has been found 
adequate for the representation of creep, stress relaxation, and strain rate 
effects encountered in a combustor liner material under service conditions 
where strain rates may vary from 10 sec to 10 sec" . However, it ap- 
pears that an exponential law is necessary if strain rates greater than about 
-2 -1 

10 sec are encountered. In particular, at higher strain rates, the power 
law expression for the inelastic strain rate and predicts values of stress in 
a constant strain rate tensile test which are too large [13, 30]. 

The model contains 14 temperature— dependent material constants. The 
model has been applied successfully for predicting creep, stress relation, 
cyclic, and thermomechanical hysteresis behavior of Hastelloy-X alloy [6] . 
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4.3 Chaboche * s Theory 


Chaboche's theory [14] is formulated on the basis of a yield surface. 
Inside the yield surface it is assumed that no inelastic deformation can take 
place. Different flow laws are thus required for loading and unloading. The 
use of a yield surface permits isotropic hardening to be modeled by an in- 
crease in the size of the yield surface rather than by an increase in the drag 
stress. Hence, in this theory, K is assumed to be constant. Initially, the 
isotropic variable Y is assumed to be zero and inelastic deformation occurs 
only when where K is a constant. As Y grows with inelastic deforma- 

tion, the yield surface expands and inelastic deformation takes place only 
when 3J>K+Y. Yield surface translation occurs as the equilibrium stress 
change s^with inelastic deformation. The growth law for the equilibrium stress 
contains a temperature rate term which allows the equilibrium stress to change 
with temperature during nonisothermal elastic excursions. The model contains 
13 constants. All of the material constants are functions of temperature and 
must be experimentally determined at each temperature of interest. 

4.4 Bodner and Partem* s Theory 

The elastic-viscoplastic theory of Bodner and Partom, (B-P) , was prob- 
ably the first "unified" set of constitutive equations without a yield cri- 
terion or loading/unloading conditions to be developed (described in the 
literature in 1968) [59]. Those equations include certain physical concepts 
provided by the work on dislocation dynamics during the 1950 's and early 
1960 * s . The equations are placed in the context of multi-dimensional con- 
tinuum mechanics which makes them capable of solving problems by analytical 
and numerical methods. One of the initial papers, in 1972, considered large 
deformations [60] and another the same year included isotropic hardening [61] 
At the present stage of development, the constitutive theory includes isotro- 
pic and directional hardening, thermal recovery of hardening, general tempera- 
ture dependence of plastic flow, and isotropic and anisotropic damage develop 
ment [7-9]. In principle, the equations could provide for the pressure de- 
pendence of plastic flow and could be expanded to include anelastic effects. 
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The B P theory differs from others of the same type in certain details, al- 
though most all other proposed theories follow the same general principles. 

The main points of difference of the B-P theory to others are the following: 

1. The basic equation for plastic strain rate as a function of stress 
and history dependent internal variables is taken to be of exponential form. 
The initial work [59], before 1970, was based on a power law relation which 
seems to have been adopted by other investigators. Properties of the exponen- 
tial function used in the B-P formulation that may be especially useful are 
its very low value, almost zero, at low stress levels and its limiting value 
of plastic strain rate. As a consequence, that particular function seems to 
be suitable to represent material behavior over a wide range of strain rates 
including the very high rates, i.e. about 10 5 sec” 1 . 

2. The scalar measure for hardening is taken to be the plastic work 
rate. However, the overall theory is not dependent on this point and other 
measures, e.g. plastic strain rate or a function of the plastic strain rate 
history, are admissable within the context of the theory. 

3. The treatment of directional hardening by Bodner and his associ- 
ates has been based on the general anisotropic form of the flow law rather 
than on the back stress" concept. An incremental isotropic form of the flow 
law has been proposed as an approximation which would simplify the numerical 
computations. In the incremental isotropic equation, the scalar coefficient 
is a function of both the isotropic and directional hardening variables which 
depend on the loading history. 

4. In the evolution equations for directional hardening, B-P theory 
uses the direct stress as the directional index of hardening. The plastic 
strain rate had been used previously for that purpose and most other theories 
still rely on that variable. 

5. Anelastic effects, i.e. reversible deformations with energy loss, 
are taken to be given by a separate parameter in the B-P theory and are not 
included in a "back stress" parameter. 

Applications of the B-P constitutive equations have been made to 
problems of steady and variable creep [57, 62], creep crack growth [63,64], 
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static and dynamic response of metal matrix composites [65,66] , dynamic frac— 
tnre mechanics [67,68], wave propagation [69], and adiabatic shear band forma- 
tion [70], The equations have been adopted in various finite element and 
finite difference computer programs, e.g. ADINA and DEPROSS [71]. 

4 .5 Stouffer and Bodner 's Theory 

This theory [18,19] is an anisotropic generalization of Bodner and 
Partom's theory for isotropic materials and is based on an anisotropic version 
of the Prandt 1— Reus s flow law. The inelastic deformation rate is related to 
the deviatoric stress by a linear transformation (a fourth order tensor) whose 
components are functions of stress, stress history, and the internal vari- 
ables. The procedure adopted by Stouffer and Bodner was to transform the 
anisotropic flow equation into six-dimensional vector form with the matrix of 
coefficients becoming a 6 x 6 second order tensor. Upon diagonal ization of 
the tensor, the strain rate component equation becomes uncoupled and allows 
the coefficients to incorporate the hardening variables in a simple manner. 

The hardening variables which represent both isotropic and directional hard- 
ening are obtained from a proposed anisotropic work hardening law written in a 
hardening/recovery format. The anisotropic formulation, however, does not 
automatically lead to plastic incompressibility and results in plastic volume 
changes [18]. The anisotropic theory was subsequently revised to enforce 
plastic incompressibility by a relatively simple recalculation of the matrix 
coefficients [19] . The model has been successfully used for the solution of a 
number of dynamic penetration and impact problems [72] • 

4 .6 Miller's Theory 

Miller's model is based on a combination of creep phenomenology and phy- 
sical mechanisms. Proposed in 1975 [10], the model uses the Garafalo hyper- 
bolic sine relation for steady state creep as a starting point and introduces 
two internal variables to describe non— steady state inelastic responses. The 
two structural variables are the drag stress and the rest stress and they re- 
present isotropic and kinematic hardening, respectively. The model has since 
been modified to include the phenomenon of solute strengthening [11] by adding 
two solute strengthening internal variables to the drag stress. The more re- 
cent version of the Miller's model also takes into account of irradiation 
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effects [73]. Excluding the solute strengthening variables, the latest ver- 
sion of Miller's model consists of four internal variables [12]. The drag 
stress and the rest stress are both decomposed into short range and long range 
components. The two isotropic hardening variables represent the hardening 
effects of obstacles such as total dislocation density, dislocation cells and 
subgrain boundaries. On the other hand, the long range and short range rest 
stress components represent the directional hardening effects due to disloca- 
tion interaction with internal stresses with long characteristic wavelength 
and short characteristic wavelength, respectively. 

The evolution rate equations of the internal variables are written in 
a hardening/ recovery format. An Arrhenius expression is used to model the 
temperature dependence of the activation energy for thermally-activated plas- 
tic flow. A dynamic recovery term is later introduced in the four-variable 
model. Both the two-variable and the four-variable Miller's model have the 
capabilities to simulate monotonic and cyclic stress-strain behavior. In both 
models, cyclic deformation behavior and the Bauschinger effect are modeled 
through the rest stress term. In the two-variable model, a linear hardening 
term is used and it results in the predictions of overly— square hysteresis 
loops. This shortcoming has been resolved by using a nonlinear hardening co- 
efficient which varies with the loading direction and becomes very large after 
a load reversal [74] • The high hardening rate in turn results in rounded and 
physical ly-real istic hysteresis loops . 

Cyclic hardening is modeled through the increase in both the drag stress 
K and the rest stress . Cyclic saturation is represented by an interaction 
between rest stress and drag stress in the hardening coefficient in the drag 
stress equation, and it requires a balance between strain hardening and dyna- 
mic recovery. The level of K at cyclic saturation depends upon the average 

absolute value of Q.. which in turn is a function of the cyclic strain ampli- 

ij 

tude . 

Steady state inelastic deformation occurs when the evolution rate equa- 
tions become zero (i.e., X. =0). This condition occurs when the hardening 

l 

rate is equal to the sum of the dynamic and thermal recovery rates. 
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The use of four internal variables in Miller's model allows the sim- 
ulations and predictions of a large variety of inelastic deformation behavior 
including complex loading histories involving strain rate changes, temperature 
transients, and hold time. The model is, however, very complex and contains a 
large number of material constants. Neglecting the constants for irradiation 
effects and the solute strengthening terms, the latest four-variable model 
contains 24 material constants. However, excluding the elastic modulus, the 
material constants are all temperature independent and need not be evaluated 
as a function of temperature. 

4.7 Kagawa and Asada' s Modification of Miller's Model 

Kagawa and Asada recently extended Miller's two-variable model to multi- 
axial form [22]. The extension is, however, quite different from Miller's 
version [73], particularly in the choice of the flow rules for treating kine- 
matic hardening. Kagawa and Asada separate the inelastic strain rate into 
isotropic and anisotropic components. The isotropic component of the inelas- 
tic strain rate is proportional to the deviatoric stress, while the anisotro- 
pic component is a linear function of the inner product of the rest stress and 
the effective stress. The evolution rate equations for the drag stress and 
the rest stress are essentially identical to Miller's model. 

The model has been used successfully to simulate and/or predict mono- 
tonic and cyclic stress-strain behavior, creep, biaxial strain ratcheting 
under cyclic loading, and equi-creep rate surface (flow potentials). 

4.8 Krieg, Swearengen, and Rohde's Theory 

This theory uses a power-law kinetic relation and is formulated in terms 
of two internal variables [13]; the drag stress is used to model isotropic 
hardening, while the back stress is used to model directional (kinematic) 
hardening. The evolution rate equations are formulated in a hardening/ re- 
covery format. However, the hardening function in the drag stress evolution 
rate equation has not been defined and is assumed to be a constant. Recovery 
of drag stress is considered a thermally-activated process and is represented 
using Friedel's climb recovery model [34]. 
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Directional (kinematic) hardening is modeled using a nonlinear hardening 
function in the back stress evolution rate equation. At large strain, the 
hardening function saturates to a constant value. The recovery of back stress 
is, again, by Friedel's climb recovery model. 

This theory does not include a dynamic recovery term. The state vari- 
ables thus saturate at large strain values when the static thermal recovery 
term is balanced by the strain hardening terms. Recovery always occurs toward 
an isotropic referenced or 'annealed" state. The recovery rate depends upon 
the magnitude by which the current state differs from the annealed state. 

The model contains nine temperature-dependent constant. Krieg et al 
[13] obtained good results between model and experiment for monotonic, cyclic, 
and creep loading conditions, but they also pointed out the inadequacy of the 
model at high strain rates. Walker [6] has evaluated and applied the model 
for predicting steady state hysteresis loops as well as thermomechanical cy- 
clic loops. 


4.9 Hart's Theory 

Hart's theory [16] employs an equilibrium stress internal variable 0^ 
with a constant drag stress E. Cyclic hardening is assumed to occur only in 

the equilibrium stress. The equilibrium stress is a function of the parameter 

♦ 

o which is referred to as the current "hardness" of the deforming materials 
by Hart [16] and his colleagues [17]. a which occurs in the evolutionary 
equations of 0^ but not in the flow law nor in the kinetic equation serves 
only to modify the equilibrium stress and may be considered as a secondary 
internal variable. 


Written in a hardening/recovery format, the equilibrium stress grows 
linearly with inelastic strain in the initial loading phase and reaches a sat- 
uration when the static thermal recovery term containing the "hardness" vari- 
able balances the linear hardening term. Since dynamic recovery terms are not 
included, the saturated value of the equilibrium stress depends on the strain 
rate. The linear hardening growth of the equilibrium stress, together with 
the rapid growth of the static thermal recovery term at large strain values. 
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produces the tril inear stress-strain response characteristic shown in in 
Delph's review paper [75]. Since the drag stress is assumed to be constant, 
the theory can model only directional (kinematic) hardening. 

Jackson et al [17] have recently modified Hart’s model by introducing 
two kinematic hardening varibles fi , , and Q . . • The hardness variable cr has 

w* ♦ 1 

also been separated into two components a 1 and The growth rates of 

and are governed by two independent evolutionary equations formulated in 

ij * * 

terms of o. and o 0 . The updated Hart model gives an improved correlation of 
the monotonic stress-strain behavior as well as cyclic hysteresis behavior at 
both large and small strain amplitudes. However, the deficiency of allowing 
cyclic hardening only in the equilibrium stress, and not in the drag stress, 
still remains. 

The introduction of an extra term in the equilibrium stress, together 
with a new hardening coefficient, o*, makes the determination of the material 
constants more difficult. The updated Hart model contains 19 material con- 
stants. The material constants which depend on temperature are not explicitly 
stated in this formulation. 

4.10 Ghosh’s Model 

Ghosh's model [21] uses only one internal state variable to describe the 
evolution of microstructures and the microscopic mechanisms of inelastic de- 
formation. The inelastic deformation processes are divided into those occur- 
ring within grains, and along grain boundaries. Thus, the total strain rate 
is the sum of the elastic strain rate, the microplastic strain rate which in- 
cludes anelasticity, the inelastic strain rate due to slip within grains and, 
finally, the inelastic strain rate due to grain boundary sliding. Different 
power laws are used to describe the kinetics of each of the inelastic strain 
rates. For deformation occurring within grains, the internal variable g des- 
cribes the internal strength within the grains. The magnitude of the internal 
strength variable g is assumed to be equal to that of the back stress. The 
evolution rate of the internal strength variable is formulated in terms of 
hardening and recovery. Both the hardening function and the recovery function 
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are dependent on the current value of the internal strength variable. Two 
modified forms of Friedel's recovery model are used to describe two stages of 
recovery which are referred to as "dynamic" and "post-dynamic" recoveries. 

The Bauschinger effect is incorporated in the model through the micro- 
oplastic strain rate term. Inelastic loading and elastic unloading are de- 
fined through the use of, and depends on whether the current stress exceeds, 
the maximum value of the internal strength variable g* obtained in the prior 
deformation histories. Nicroplasticity occurs when |ol<g* and is modeled 
through the back stress and the microplastic strain, which is allowed to be 
stored 1 and "recovered" in a time-dependent manner. The inelastic strain 
rate due to grain boundary sliding is taken as a Newtonian viscous type of 
flow. 

The model contains 13 constants and they must be evaluated as a function 
of temperature. The model has been used to simulate monotonic and cyclic 
stress-strain behavior, creep, and stress relaxation. The model has been for- 
mulated and applied in one-dimensional form only. Extension of this theory to 
multidimensional form appears to be difficult because of the need to keep 
tract of g* to define loading and unloading conditions. 
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5.0 PREDICTIVE AND SIMULATIVE CAPABILITIES OF UNIFIED 
CONSTITUTIVE THEORIES 

The unique characteristic that distinguishes the unified theories 
from constitutive theories based on classical plasticity or creep approaches 
is the ability of the unified theories to predict or simulate material re- 
sponses under monotonic, cyclic, creep, and stress relaxation loading con- 
ditions by using the same set of internal variables and material constants. 

At present, not all unified theories can perform all the four categories of 
inelastic behavior. The capacities of some of the models have been demon- 
strated through numerical exercises, but they have not been severely tested 
by comparing with experimental data. Four of the more established models 
which have been successfully applied for simulating and/or predicting mono- 
tonic, cyclic, creep, and stress relaxation behavior are those of Robinson 
[15], Walker [6], Bodner-Partom [7-9], and Miller [10-12]. Robinson's model 
is based on a yield condition and utilizes loading and unloading criteria, 
while the latter three do not. The kinetic equations commonly used in unified 
theories without a yield surface or flow potential are based on the power- law 
exponential, and hyperbolic sine functions; these kinetic equations are repre 
sented in Walker, Bodner-Partom and Miller theories, respectively. The simu- 
lative and predictive capabilities of these four unified theories are illus- 
trated below. 

5 .1 Monotonic Stress— Strain Behavior 

All unified theories are capable of reproducing the monotonic stress- 
strain curve. Most unified theories use the monotonic stress-strain data as 
part of the data base from which model constants are evaluated. Figure 9 
shows an experimental uniaxial tensile stress-strain curve of Hastelloy-X 
deformed at a strain rate of 1.3 x 10 sec at 922 K and model simulation 
using the Bodner-Partom theory. The computed curve includes contributions 
from both work hardening and thermal recovery. The use of the exponential 
function in the kinetic relation allows Bodner-Partom' s theory to simulate 
stress-strain response over a wide range of strain rates including those in 
the dynamic range. Figure 10 shows the calculated and experimental static, 
dynamic, and incremental shear stress-strain curves for copper at 298°F [54] 
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It can be seen in Figure 10 that Bodner— Partom' s model duplicates well the 

“4 — i 7 _i 

stress-strain of copper at both low (2 x 10 sec ) and high (3 x 10 sec ) 
strain rates. 

5 .2 Cyclic Stress-Strain (Hysteresis) Behavior 

The Bauschinger effect is represented in most unified theories by a 
kinematic or directional hardening internal variable. Cyclic hardening, how- 
ever, can be represented by increases in the isotropic hardening variable, the 
kinematic hardening variable, or both. These different types of cyclic hard- 
ening behavior are illustrated in Figure 11 for Walker's theory, which uses an 
equilibrium stress for modeling kinematic hardening, and in Figure 12 for 
Bodner— Partom ' s theory which does not use an equilibrium stress. Cyclic hard- 
ening depicted in the hysteresis loops of Figures 11a and b are due to the 
increases of the equilibrium stress and the drag stress, respectively. On the 
other hand, cyclic hardening in Bodner— Partom theory are the consequences of 
the increases in (a) the directional component, and (b) both the isotropic and 
directional components of the internal variable Z, as shown in Figure 12a and 
b, respectively. Despite the different approaches in treating directional 
hardening, both the Walker and Bodner— Partom theories yield realistically 
rounded hysteresis loops. This is demonstrated in Figure 13 which shows com- 
parisons of theoretical calculations of saturated hysteresis loops based on 
these two theories and experimental data of Hastelloy-X deformed at five dif- 
ferent strain rates at 1144°K. 

The use of different growth laws of the equilibrium stress in different 
regions of stress space allows Robinson's model to reproduce rounded hyster- 
esis loops. Examples of cyclic saturated hysteresis loops calculated using 
Robinson's model are compared with experimental results of 2-l/2Cr-lMo steel 
in Figure 14 [55] . 

5 .3 Creep Responses 

Most of the unified models can predict or simulate primary and secondary 
creep responses of material subjected to constant load or stress. Steady 
state creep rates are predicted by these unified models to occur when the evo- 
lutionary rates of the isotropic and/or directional hardening variable vanish 
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(a) (b) 


FIGURE 11. HARDENING OF HYSTERESIS LOOP IN WALKER'S MODEL. 

(a) Hardening due to equilibrium stress; (b) 
Hardening due to drag stress, from [6], 




FIGURE 12. HARDENING OF HYSTERESIS LOOP IN BODNER-PARTOM'S MODEL. 

(a) Hardening due to Z&; (b) Hardening due to ~L\ and 

ZD. 
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FIGURE 13. COMPARISON OF EXPERIMENTAL AND FUNCTIONAL THEORY PREDICTIONS FOR HASTELLOY-X AT 871 °C (1600°F) 
(a) Experimental data, (b) Walker's Theory [6], and (c) Bodner-Partom's Theory. 
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as the hardening terms are balanced by the thermal recovery. Examples of cal- 
culated steady state creep rates under constant stress and comparison with 
experimental data are shown in Figure 15 and Figure 16 for Robinson's model 
[53] and Bodner-Partom' s model, respectively. According to the unified theo- 
ries, the steady state creep rate is a function of stress and temperature 
only; it should not depend on the loading histories. This is demonstrated in 
both experimental data and predictions by Miller's model in Figure 17 which 
shows the response of a creep specimen subjected to a sudden decrease in ap- 
plied stress from 27.6 MPa to 13.8 MPa after a creep strain of 16 percent. 
Miller's calculation [56] indicates that the 16 percent strain at 27.6 MPa 
results in an internal variable which is larger than the steady-state one at 
13.8 MPa; hence, when the stress is decreased, the creep rate first drops 
instantaneously but then gradually increases as recovery reduces the value of 
the internal variable to one which is characteristic of the lower stress. 

The instantaneous creep response of material subjected to an arbitrary 
loading history depends on the current values and the growth laws of the in- 
ternal variables. For unified theories based on an equilibrium stress and 
without a flow potential, the creep response would, according to the flow law, 
depend on the difference between the current stress and the equilibrium 
stress. Different creep responses would result from the same imposed stress 
(Points A and D in Figure 18) on the loading and unloading branches of a cy- 
clically-saturated hysteresis loop because of differences in the equilibrium 
stresses, as demonstrated in Figure 18 which shows both Walker's model pre- 
diction and experimental data of Hastelloy-X at 1144°K [6]. Walker T s model 
also predicts that a compressive creep strain can occur under a constant ten- 
sile stress if the equilibrium stress is algebraically larger; this reverse 
creep behavior is characteristic of unified theories which are based on the 
equilibrium stress and without a yield surface. 

5.4. Stress Relaxation Response 

The behavior of unified constitutive models under stress relaxation is 
analogous to the creep behavior. Under a constant strain condition, the re- 
laxation rate would, again, depend on the current values of the internal vari- 
ables and on the growth laws which describe their changes with time and in- 
elastic deformation. Stress relaxation calculations based on Bodner-Partom' s 
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FIGURE 16. STEADY CREEP RATES AS A FIGURE 17. MILLER'S MODEL PREDICTION 

FUNCTION OF STRESS SIMU- COMPARED WITH EXPERIMENTAL 

LATED BY BODNER-PARTOM'S DATA FOR A CREEP TEST WITH 

MODEL A SUDDEN DECREASE IN AP- 

PLIED STRESS, FROM [28] 
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x 10-3 SEC-1 WITH A STRAIN AMPLITUDE OF LEVELS. The calculated 

± .4%. Calculated curves are generated curves are generated with 

with Walker's theory, from [6]. Bodner-Partom' s theory, 

from [57]. 



model and Walker's theory are compared with experimental data of Rene 95 [57] 
and Hastelloy-X [6] in Figure 19 and Figure 20, respectively • Figure 20 shows 
that cyclically saturated Hastelloy-X subjected to compressive strain hold at 
1144° X relaxes from an initial compressive stress to a tensile one after ap- 
proximately 3 seconds. This behavior is predicted qualitatively though not 
quantitatively by Walker's model on the basis of a positive equilibrium 
stress, as shown in Figure 20. 

5 .5 Thermomechanical Response 

The behavior of unified constitutive theories under thermomechanical 
cycling depends critically on the change of material constants with tempera- 
ture. In particular, the shape of the predicted thermomechanical loop is 
sensitive to the growth of the kinematic hardening variable (the equilibrium 
stress) with temperature. Walker has also found that it is necessary to 
include a temperature rate term in the evolutionary equation such that the 
equilibrium stress can change during nonisothermal "elastic" excursions. 
Walker's model prediction of thermomechanical loop of Hastelloy-X is shown in 
Figure 21 [6]. 

5 .6 Multiaxial Behavior 

All the unified theories utilize single-valued kinetic equations formu- 

2 ' 2 

lated in terms of either SJ^/X or 3J^/X . For a constant value of the inter- 
nal X, these kinetic equations predict a locus of constant inelastic strain 
rate invariant = constant) in stress space; the shape of the predicted 

"yield surface" or "flow potential" is identical to von Mises yield function 
and is described by 

3J 2 - C 2 * 2 =0 or 3J 2 ’ - C 2 K 2 = 0 

where C is a constant. This is illustrated in Figure 22 which shows that for 

proportional loading the "yield surface" predicted by Bodner-Partom' s model 

is identical to the von Mises one, and the inelastic strain rate vectors are 

all normal to the "yield surface." For unified models formulated based on the 

* 

equilibrium stress. Equation (10) remains valid. With J 2 given by Equation 
(8b), the size of the "yield surface" is given by CK, while the center of the 
"yield surface" is at and translates according to the evolution rate of 
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. This is demonstrated in Figure 23 which shows the equi-creep rate 
surfaces calculated for proportional loading under combined tension and tor- 
sion using a modified Miller theory [22]. The broken line shows the initial 
creep surface, and the solid line depicts the subsequent surface after a 
3 x lO” 1 prestrain in tension. The equivalent creep rate for both surfaces is 
1 x 10 ^ sec The subsequent creep surface shows both an expansion and a 
translation of the origin in the prestrain direction* 
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6.0 SUMMARY AND CONCLUSIONS 


1. A review of more than ten time-dependent elastic-viscoplastic 
constitutive theories indicates that these theories differ in 
the choice of flow law, kinetic equation, and evolutionary 
equation of the internal variables. 

2. The unified approach treats all aspects of inelastic deformation 
including plasticity, creep, and stress relaxation using the same 
set of flow law, kinetic equation, and internal variables. 

3. The unified constitutive theories satisfy the uniqueness and 
stability criteria imposed by Drucker's postulate for rate in- 
dependent stable plastic flow and Ponter's inequalities for con- 
stitutive theories based on internal variables. 

4. The unified theories can be formulated either with or without the 
use of a yield criterion. Three basic flow laws are identified in 
theories without a yield criterion. For theories with a yield 
criterion, the associated flow law is derived from the yield 
function or the flow potential. 

5. Three different formulations of the kinetic equations are identified, 
and they include the exponential, power law, and hyperbolic sine 
functions. The exponential formulation gives a limiting inelastic 
strain rate and appears to give better results for high strain rate 
applications. 

2 

6. All forms of kinetic equations reviewed are functions of 3J 0 /K (or 

r 2 2 
3J jl^L ) and result in ’yield surfaces’’ and equi-creep rate sur- 
faces which are described by the von Mises criterion. 

7. The number of internal variables used in the unified theories range 
from one to six. Most unified theories use two internal variables, 
one to represent isotropic hardening and one to present kinematic 
or directional hardening. The measure of hardening is either the 
inelastic strain rate or the inelastic work rate. 
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8. Directional (kinematic) hardening can be modeled with or without the 
use of an equilibrium stress. The directional index of kinematic 
hardening can be based on either the inelastic strain rate or the 
direct stress. 

9. Material constants in the unified models are necessarily temperature- 
dependent and required to be evaluated at the temperatures of inter- 
est. The change of material constants with temperature has a drastic 
effect on the shape of the thermomechanical hysteresis loops. There 
are indications that a temperature rate term is also required in the 
unified theories. 

10. All of the unified theories which are reviewed do not automatically 
predict additional cyclic hardening under nonproportional loading 
paths. Additional terms are needed in the unified theories to in- 
clude such hardening behavior. 

11. The equilibrium-stress-based unified theories can describe reverse 
creep and/or reverse stress relaxation behavior without further 
modifications. Unified models which are not based on the equilibrium 
stress would require modification by adding an anelastic term in order 
to take into account these types of behavior. 

12. The unified constitutive equations are stiff but can be integrated using 
either explicit or implicit methods. The choice of a particular inte- 
gration scheme in finite element analysis appears to be code and problem 
dependent . 
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